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ABSTRACT 


The  threat  of  collision  between  an  asteroid  or  a  comet  and  the  Earth  has  been  well 
docxunented.  Mitigation  of  such  a  threat  can  be  accomplished  by  destruction  of  the  threat 
or  by  perturbing  the  threat  object  into  a  safe  orbit  Following  a  summary  of  proposed 
mitigation  techniques,  this  thesis  investigates  the  impulse  required  to  safely  perturb  a 
threatening  Earth-Crossing  Asteroid  (ECA).  While  previously  published  analysis 
included  only  two-body  approximations  to  the  impact  geometry,  this  thesis  adds  the 
effect  of  the  Earth's  gravitational  field  to  more  closely  approximate  reality.  The  results 
indicate  that  third-body  effects  are  strongest  on  ECA's  in  a  nearly  circular  heliocentric 
orbit,  where  the  minimum  required  AV  can  be  several  times  larger  than  that  calculated 
using  two-body  approximations.  To  determine  the  minimum  AV  required  for  mitigation, 
MATLAB®  's  sequential  quadratic  programming  (SQP)  algorithm  is  applied  to  a 
constrained  optimization  problem.  Third-body  effects  were  added  to  a  previously 
published  two-body  optimization  by  modifying  the  boundary  conditions.  With 
knowledge  of  the  minimum  AV  requirements,  the  capability  of  current  impulsive 
mitigation  technology  is  analyzed.  For  asteroids  of  median  density  in  co-planar  orbits,  a 
single  24  Mt  nuclear  explosive  impulse  applied  earlier  than  3  years  before  impact  can 
effectively  mitigate  a  threat  with  a  diameter  of  6  km.  The  capability  significantly 
decreases  with  shorter  warning  times. 
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INTRODUCTION 


Who  knows  whether,  when  a  comet  shall  approach  this  globe  to  destroy  it, 
as  it  often  has  been  and  will  be  destroyed,  men  will  not  tear  rocks  from 
their  foundations  by  means  of  steam,  and  hurl  mountains,  as  the  giants  are 
said  to  have  done,  against  the  flaming  mass?  -  And  then  we  shall  have 
traditions  of  Titans  again,  and  of  wars  with  Heaven. 

-  Lord  Byron,  1822  (Medwin,  1824,  p.  185j 

Although  steam  is  no  longer  the  motive  force  of  choice,  the  defense  that  Lord 
Byron  dreamed  of  150  years  ago  is  finally  possible  today.  For  the  first  time  in  the 
history  of  humanity,  the  technology  exists  which  could  potentially  protect  the  Earth 
against  an  impact  fi*om  an  asteroid  or  comet.  In  addition,  the  threat  of  impact  from  a 
cosmic  body  that  Lord  Byron  addressed  has,  within  the  last  20  years,  become  widely 
accepted.  However,  while  the  defensive  technology  exists,  an  active  planetary  defense 
plan  does  not.  At  this  incipient  stage  in  researching  impact  hazards,  there  are  still  many 
challenges  to  be  met. 

Before  considering  the  technical  challenges  to  a  defense  plan,  one  must  first 
answer  the  question:  what  is  the  threat?  Cosmic  bodies  in  orbits  that  bring  the 
possibility  of  impact  with  the  Earth  are  commonly  termed  Near  Earth  Objects  (NEOs). 
Generally,  NEOs  are  comets  or  asteroids  that  have  a  perihelion  radius  of  less  than  1.3 
AU.  These  are  orbits  that  have  the  potential  to  be  perturbed  by  the  gravitational  force  of 
the  Earth  and  other  planets  into  orbits  that  impact  the  Earth  (Cheng,  1994,  p.  651). 

Lord  Byron's  quotation  refers  to  comets,  which  are  distinguished  from  asteroids 
by  outgassing  activity  or  the  existence  of  a  tail  (Cheng,  1994,  p.  651).  Comets  are 
believed  to  be  composed  of  ice  and  rock  (Elder,  1997,  p.  11),  and  can  approach  in  such 
highly  eccentric  orbits  that  they  can  impact  the  Earth  at  velocities  greater  than  50  km/s 
(Morrison,  1994,  p.  61).  The  composition  of  asteroids  varies  widely  fi-om  metal  to  the 
same  ice  and  rock  foimd  in  comets  (Morrison,  1994,  p.  61).  Because  of  the  lower 
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eccentricity  of  their  orbits,  asteroids  generally  impact  the  Earth  at  velocities  on  the  order 
of  20  km/s  (Morrison,  1994,  p.  61). 

The  technical  challenges  to  a  NEO  defense  plan  can  be  grouped  into  two 
categories:  hazard  detection  and  hazard  mitigation. 

Current  theory  suggests  that  only  a  small  percentage  of  all  potentially  hazardous 
Near  Earth  Objects  (NEOs)  have  been  discovered.  It  is  believed  that  all  Earth-Crossing 
NEOs  with  diameters  of  6-12  km  have  been  discovered.  With  smaller  diameters,  the 
detection  completeness  decreases  dramatically:  about  35%  with  diameters  of  3-6  km,  15% 
with  diameters  from  2-4  km,  and  7  %  with  diameters  from  1-2  km.  This  low  detection 
completeness  implies  that  a  significant  number  of  objects  with  the  potential  of  causing  a 
global  disaster  have  not  yet  been  detected.  Some  specific  challenges  for  detection  are: 
correcting  a  bias  in  discoveries  that  results  from  the  majority  of  searches  being  conducted 
in  the  Earth's  northern  hemisphere,  improving  discovery  follow-up  to  achieve  a  reliable 
orbit  calcvdation  (Carusi,  1994,  p.  145),  decreasing  the  minimum  possible  detection 
magnitude,  and  coherently  processing  detection  data  (Bo well,  1994,  p.  194). 

Once  detected,  the  orbital  parameters  of  the  NEO  must  be  accurately  determined. 
Errors  in  the  estimation  of  a  NEO's  orbital  parameters  can  propagate  over  time  to  produce 
a  drastically  inaccurate  estimate  of  impact  probability  or  time.  A  thorough  imderstanding 
of  this  error  propagation  is  critical  to  the  planning  of  a  defense  mission.  Another 
challenge  is  to  model  the  non-gravitational  effects  in  comets  caused  by  rocket-like 
outgassing  (Yoemens,  1994,  p.  257)  so  that  accurate  impact  probabilities  can  be 
calculated. 

This  thesis  concentrates  on  some  of  the  challenges  to  the  mitigation  of  a  NEO 
hazard.  Chapter  II  presents  a  detailed  study  of  mitigation  ideas.  The  thesis  then 
continues  with  the  goal  of  answering  two  questions: 

1 .  What  is  the  minimum  AV  required  to  deflect  a  threat-NEO  into  a  safe  orbit? 
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2.  What  is  the  maximum  NEO  size  that  can  be  deflected  with  a  single  defense 
mission  using  current  technology? 


The  challenges  addressed  in  this  thesis  are  only  a  fraction  of  those  that  must  be 
overcome  before  a  Planetary  Defense  system  can  be  fielded. 
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IL  STATEMENT  OF  THE  PROBLEM 


A  HAZARD  MITIGATION  TECHNIQUES 

All  proposed  NEO  hazard  mitigation  techniques  (Simonenko,  1994;  Wood,  1994; 
Ahrens,  1994;  Solem,  1994;  Canavan,  1994;  Meissinger,  1995;  Melosh,  1994;  Shafer, 
1994)  attempt  to  accomplish  one  of  two  goals:  either  destroy  the  NEO,  producing 
fragments  that  will  not  inflict  damage  if  they  impact  the  Earth,  or  deflect  the  NEO  into  an 
orbit  that  is  no  longer  hazardous.  Both  of  these  techniques  require  a  controlled  coupling 
of  energy  to  the  NEO.  Much  of  the  analysis  of  hazard  mitigation  techniques  thus  far  has 
revolved  around  determining  the  best  source  of  energy  for  hazard  mitigation  (Canavan, 
1994;  Ahrens,  1994;  Simonenko,  1994;  Shafer,  1994;  Solem,  1994;  Willoughby,  1994; 
Melosh,  1994). 

As  with  all  space  missions,  the  mass  of  the  payload  is  a  critical  driver  in  the 
selection  of  energy  source.  Table  2.1  gives  a  comparison  of  the  specific  energies  of  many 
potential  energy  sources.  If  energy  source  selection  were  based  solely  on  payload  mass, 
nuclear  explosive  would  obviously  be  the  best  choice. 


High  Explosive 

6MJ/kg 

Kinetic  Energy  (10  km/s) 

50  MJ/kg 

Nuclear  Explosive 

4*10^MJ/kg 

Table  2.1  Specific  Energies  (after  Shafer,  1994) 


It  is  instructive  at  this  point  to  discuss  terminology  used  in  describing  nuclear 
explosives.  In  most  technical  applications,  the  size  of  a  nuclear  explosive  is  expressed  as 
the  amount  of  energy  the  explosive  can  produce.  This  energy  is  usually  measured  in 
terms  of  an  equivalent  mass  of  TNT  that  would  produce  the  same  explosive  power. 

Using  this  convention,  then,  1  MT  =  4.2  x  10*^  Joules  (Morrison,  1994,  p.  61).  Table  2.2 
presents  the  mass  of  a  nuclear  explosive  charge  based  on  its  yield. 
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Yield 

Mass 

1  Mt 

0.5  ton 

10  Mt 

3  to  4  ton 

100  Mt 

20  to  25  ton 

Table  2.2  Yield  Versus  Mass  for  Nuclear  Explosive  Devices  (after  Simonenko,  1994,  p. 

931) 

However,  mass  is  only  one  of  several  factors  that  must  be  considered  in  choosing 
an  energy  source.  The  following  is  a  compilation  of  proposed  NEO  hazard  mitigation 
techniques,  based  on  the  mitigation  goal. 

1.  Fragmentation  and  Dispersal 

If  sufficient  energy  is  delivered  to  the  NEO  to  fragment  it,  the  hazard  could  be 
mitigated  in  several  ways.  First,  differential  velocities  in  the  individual  fragments  of  the 
destroyed  NEO  would  result  in  a  dispersed  debris  cloud.  Some  small  scale  fragmentation 
experiments  suggest  that  if  a  NEO  were  destroyed  1  orbit  prior  to  impact,  as  little  as  1% 
of  the  original  NEO  mass  would  impact  the  Earth  (Ahrens,  1994,  p.  919). 

The  second  mitigation  effect  of  fragmenting  the  NEO  is  to  reduce  the  size  of  the 
meteor  that  reaches  the  Earth.  In  order  to  be  effective,  a  mitigation  effort  would  be 
required  to  reduce  the  impactor  size  to  the  point  that  it  Would  no  longer  generate  damage 
on  impact  with  the  Earth.  The  atmosphere  can  provide  a  great  deal  of  protection  from 
some  meteors,  but  its  effects  are  highly  dependent  on  the  impact  velocity  and  the 
composition  of  the  NEO.  Figure  2.1  shows  the  relation  between  the  composition  of  a 
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NEO,  its  impact  velocity,  its  size,  and  the  amount  of  destruction  that  can  result  following 


atmospheric  ablation  and  impact. 


Figure  2.1  The  radius  of  destruction  around  the  impact  point  due  to  the  atmospheric  blast 

wave  (from  Hills,  1993,  p.  1 133). 


But  what  is  the  best  way  to  fragment  the  NEO?  For  the  mass  reasons  discussed 
above,  nuclear  explosives  are  a  very  attractive  option  for  the  fragmentation  energy  source. 
Smaller,  more  loosely  bound  NEOs  could  be  destroyed  by  a  surface  explosive.  However, 
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larger  and  more  solid  NEOs  would  require  that  the  explosive  be  buried.  Burial  would 
provide  significantly  greater  coupling  of  energy  to  the  NEO  (Ahrens,  1994,  p.  917). 

The  requirement  for  charge  burial  adds  some  complexity  to  the  mitigation  effort. 
Wood  (1994)  suggests  some  ways  to  effect  this  burial.  A  penetration  depth  of 
approximately  10  meters  could  by  accomplished  by  use  of  hyper- velocity  (kinetic 
energy)  blasting  prior  to  arrival  of  the  nuclear  charge.  For  larger  NEOs,  a  string  of  20-100 
0.5  MT  nuclear  explosives  could  be  used  for  penetration.  Once  the  desired  depth  is 
achieved,  a  20-300  MT  charge  can  be  directed  into  the  void  and  detonated  to  achieve  final 
NEO  destruction. 

Ahrens  and  Harris  (1994,  p.  922)  estimate  that  a  charge  buried  at  optimal  depth 
can  destroy  a  hazardous  NEO  as  summarized  in  Table  2.3: 


NEO  Diameter 

Required  Nuclear  Explosive  Yield 

0.1  km 

800  kg 

1  km 

22  Kt 

10km 

0.6  Gt 

Table  2.3  Nuclear  Explosive  Yield  required  to  fragment  hazardous  NEOs  (after  Ahrens, 

1994,  p.  922) 

While  analysis  indicates  that  nuclear  explosives  can  successfully  fragment  a  NEO, 
there  are  hazards  associated  with  this  energy  source.  First  is  the  concern,  presented  by 
Meissinger  (1995,  p.  16),  that  the  NEO  fragments  which  reach  ftie  Earth  following  nuclear 
destruction  could  be  radioactive.  The  most  significant  problem  with  using  nuclear 
weapons  in  NEO  mitigation  is  the  tremendous  political,  sociological,  economic,  and 
management  aspects  of  the  weapons. 

The  alternatives  for  NEO  fragmentation  use  the  kinetic  energy  inherent  in  the 
closing  speeds  between  the  NEO  and  an  interceptor.  Wood  (1994)  proposes  a  "Jack- 
Hammer"  rubblization  technique  in  which  waves  of  hyper-velocity  projectiles  are  directed 
towards  the  threat  object.  Each  projectile  vaporizes  a  small  portion  of  the  object,  leaving 
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a  hole  for  the  next  projectile  to  enter.  The  cumulative  effect  of  the  impact  of  many 
projectiles  is  several  deep  holes  and  a  somewhat  fragmented  NEO. 

Woods'  team  further  refines  this  concept  (Teller,  1995)  into  Hypervelocity 
Projectile  Array-Sheets  and  Hypervelocity  Projectile  Lattices.  The  Sheets  are  a  two 
dimensional  plane  with  multiple  projectiles  connected  to  each  other.  The  connection 
helps  to  ensure  that  follow-on  projectiles  will  be  directed  to  impact  spots  of  previous 
projectiles.  The  Lattice  combines  several  Sheets  into  a  three  dimensional  pulverization 
effort. 

Thus,  there  are  several  options  for  fragmentation  and  dispersal  of  a  threat-NEO. 
There  are,  however,  arguments  against  using  fragmentation  as  a  mitigation  goal.  In  order 
for  fragmentation  to  significantly  reduce  the  amount  of  NEO  mass  that  eventually 
impacts  the  Earth,  the  NEO  must  be  destroyed  early.  Thus,  Ahrens  and  Harris  conclude 
that  "...fragmentation  is  likely  to  be  a  safe  choice  only  for  long  lead-time  response 
(decades)  or  for  relatively  small  bodies  where  the  fragments  may  be  allowed  to  hit  the 
Earth  (Ahrens,  1994,  p.  919)."  If  the  target  were  pulverized  late,  it  may  make  the 
problem  worse.  The  resulting  fragments  would  present  multiple  targets  to  any  follow-on 
mitigation  effort  (Solem,  1994,  p.  1028). 

■  2.  Orbital  Deflection 

Many  NEO  orbital  deflection  schemes  work  by  accelerating  mass  away  from  the 
NEO,  thus  changing  its  momentum.  Chemical  rockets,  which  have  seen  extensive  use  in 
space  for  changing  the  momentum  of  man-made  satellites,  are  only  appropriate  for  NEOs 
with  diameters  of  less  than  100  meters  (Wood,  1994,  p.  1 1).  The  mass  of  fuel  required 
to  rendezvous  (by  matching  velocity)  with  the  threat  and  then  to  displace  it  would  be 
otherwise  be  excessive. 

Another  option  is  to  use  a  mass  driver.  A  mass  driver  uses  electromagnetic  forces 
to  accelerate  buckets  containing  material  from  the  surface  of  the  NEO  (Melosh,  1994,  p. 
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1117).  Thus,  a  mass  driver  has  an  advantage  over  chemical  rockets  in  that  it  avoids  the 
need  to  transfer  propellant  mass  from  the  Earth.  While  the  technology  for  mass  drivers 
has  been  available  since  the  mid-1970's  (Melosh,  1994,  p.  1 1 17),  they  are  nonetheless 
very  complex.  Meissinger  (1995)  points  out  that  mass  drivers  would  require  large-scale, 
complex  robotic  operations  for  construction  and  processing  of  the  NEC's  soil.  Also,  the 
mass  driver  concept  is  critically  dependent  on  the  soil  conditions  of  the  NEC  (Meissinger 
1995).  At  this  point,  little  is  known  about  the  composition  or  soils  characteristics  of 
asteroids  and  comets. 

Just  as  the  kinetic  energy  of  the  threat-NEO  was  used  in  NEC  destruction 
schemes,  there  are  schemes  to  use  this  kinetic  energy  for  NEC  deflection.  Interestingly, 
for  kinetic  energy  mitigation  schemes,  Canavan  (Canavan,  1994,  p.  102)  suggests  that 

"...the  mass  required  scales  on  -r,  so  faster  NEOs  present  less  of  a  threat  because  of  their 

V 

higher  specific  energy."  In  these  schemes,  a  projectile  impacts  the  threat-NEO  and  the 
crater  ejecta  combines  with  the  momentum  of  the  projectile  to  produce  an  impulse. 
Analysis  conducted  by  Ahrens  and  Harris  (1994,  p.  904)  suggests  that  the  energy  from  a 
kinetic  impact  is  coupled  much  more  efficiently  than  the  energy  from  nuclear  weapons 
detonated  at  the  surface  of  the  NEC.  Figure  2.2  shows  an  estimate  of  the  capability  of 
kinetic  energy  deflectors.  The  three  lines  show  the  capability  of  an  impactor  with  a 
diameter  of  Im,  10m,  and  100m.  Melosh  (Melosh,  1994,  p.  1115)  even  suggests  a 
"billiards  shot"  scenario,  in  which  a  small  asteroid  is  displaced  into  a  larger  asteroid. 

Wood  (1994)  also  refines  the  kinetic  energy  approach  into  a  "hypervelocity  sandrblaster." 
In  this  scheme,  a  steady  stream  of  projectiles  are  directed  at  the  threat-NEO,  combining 
their  effect  to  produce  the  required  momentum  change. 
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Figure  2.2  Capability  of  Kinetic  Energy  Deflectors  (from  Melosh,  1994,  p.  1116) 

The  next  potential  technique  for  orbital  deflection  is  the  use  of  nuclear  explosives. 
Momentum  change  can  be  imparted  on  a  threat-NEO  using  a  nuclear  explosive  in  one  of 
two  fundamental  ways:  a  stand-off  radiative  explosion  or  a  siuTace  explosion.  For  the 
stand-off  case,  the  explosive  is  designed  to  have  a  substantial  fraction  of  its  yield  as 
energetic  neutrons  and  gamma  rays.  These  rays  irradiate  the  surface  of  the  NEO  that  is 
exposed  to  the  explosion.  The  irradiated  surface  then  expands  and  ablates  away  from  the 
NEO.  As  the  material  ablates  away,  the  momentum  of  the  NEO  changes.  The  surface 
explosion  works  in  much  the  same  way  as  a  kinetic  defense  scheme  by  creating  a  crater. 
The  ejecta  departing  from  this  crater  creates  a  momentum  change  in  the  NEO.  (Ahrens, 
1994,  p.  910). 

Not  surprisingly,  it  is  not  yet  clear  which  nuclear  explosive  method  has  the 
greatest  mitigation  capacity.  Surface  nuclear  explosions  are  generally  considered  more 
effective  in  coupling  energy  into  the  NEO.  Canavan  (Canavan,  1994,  p.  105)  concludes 
that  "...the  energy  required  for  stand-off  is  greater  than  that  for  slightly  subsurface  bursts 
by  a  factor  of  about  40."  Solem  (1994,  p.  1032)  reaches  the  same  conclusion,  stating  that 
a  stand-off  deflection  would  require  a  substantial  increase  in  interceptor  mass  over  a 
surface  burst.  Ahrens  and  Harris  (1994,  p.  917),  in  contrast,  propose  that  "...surface 
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explosions  appear  to  be  not  substantially  better  than  radiative  stand-off  explosions,  in 
deflecting  NEOs."  One  significant  disadvantage  to  the  surface  explosion  is  the  risk  of 
inadvertently  fragmenting  the  threat  object.  As  pointed  out  previously,  a  fragmented 
NEO  presents  multiple  targets  for  any  follow-on  mitigation  effort.  But  there  is  fairly 
consistent  agreement  that  the  interceptor  weight  for  a  nuclear  mitigation  effort  is  several 
orders  of  magnitude  less  than  that  of  a  kinetic  energy  effort  (Solem,  1994,  p.  1032). 

Lasers  also  have  a  potential  use  in  NEO  hazard  defense.  With  this  technique,  a 
high  energy  laser  is  directed  at  the  surface  of  the  threat-NEO.  This  creates  a  thermal  flux 
which  ablates  the  surface,  much  like  the  radiative  heating  does  in  stand-off  nuclear 
explosions.  One  estimate  suggests  that  a  laser  output  of  1  GJ/s  for  12  uninterrupted  days 
could  match  the  energy  fluence  of  a  1  Mt  nuclear  burst.  (Shafer,  1994,  p.  965) 

Solar  sails  have  also  been  proposed.  Solar  sails  use  solar  radiation  pressure  to 
provide  a  motive  force.  This  force  is  small,  but  consistent  over  long  periods  of  time.  The 
long  build  up  could  provide  a  significant  change  in  momentum  of  a  hazardous  NEO. 
However,  Melosh  points  out  that  "...truly  enormous  structures  are  necessary  to  deflect 
asteroids  in  the  1  to  10  km  diameter  range."  Figure  2.3  demonstrates  the  required  size. 
The  lines  are  for  solar  sails  of  10  km,  100  km,  and  1000  km  diameters.  The  technology 
does  not  yet  exist  to  build  such  enormous  solar  sails  at  great  distances  from  the  Earth. 
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Advance  time,  yr 

Figure  2.3  Capability  of  Solar  Sails  (from  Melosh,  1994,  p.  1 120) 

One  of  the  more  interesting  non-nuclear  approaches  is  a  solar  collector.  Melosh 
(Melosh,  1994,  p.  1 120)  has  done  some  extensive  analysis  of  this  approach  and  concludes 
that,  for  non-nuclear  alternatives,  it  is  "...an  approach  that  is  arguably  better  than  any 
other  previously  proposed."  This  scheme  also  uses  a  solar  sail,  but  instead  of  using  the 
solar  sail  directly  as  a  motive  force  the  sail  is  used  to  focus  sunlight  onto  the  surface  of 
the  NEO.  The  surface  is  thus  vaporized,  and  the  ablation  of  surface  material  produces  a 
momentum  change.  Melosh  indicates  that  a  1  km  diameter  solar  collector,  which  can  be 
launched  by  the  space  shuttle,  can  deflect  asteroids  up  to  3.4  km  in  diameter  after 
operating  for  a  year  (Melosh,  1994,  p.  1125). 

B.  PREVIOUS  ORBITAL  DEFLECTION  ANALYSIS 

The  frmdamental  goal  of  this  thesis  is  to  investigate  the  effects  of  third-body 
gravitation  on  orbital  deflection  requirements.  Previous  research  has  simplified  die 
problem  by  assuming  two-body  orbital  mechanics  between  the  Sun  and  the  threat-NEO. 
This  assumption  neglects  perturbations  due  to  the  Earth's  gravity.  While  these 
perturbations  may  not  be  present  until  the  terminal  phase  on  the  impact  scenario,  they 
affect  both  long  and  short  warning  time  analyses  as  described  below. 
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1.  Long  warning  times 


Most  analyses  of  long  warning  time  threats  have  concentrated  on  changing  the 
phasing  of  the  NEO.  Ahrens  (1994,  p.  903)  approximated  the  required  AV  for  deflection 
by  comparing  the  position  of  the  perturbed  NEO  orbit  with  that  of  the  unperturbed  orbit. 
He  then  calculates  the  velocity  change  required  to  produce  a  displacement  of  1  Earth 
radius  (R^.  Investigation  of  third  body  effects  and  patched  conic  theory  indicates  that  a 
displacement  of  1  R®  may  be  insufficient. 

Earlier  work  done  at  the  Naval  Postgraduate  School  (Knudson,  1995;  Park,  1997; 
Elder,  1997)  refines  the  estimates  of  the  required  AV.  Ail  of  these  analyses,  however, 
begin  with  a  two-body  approximation.  One  would  like  to  know  how  much  third-body 
gravitation  will  affect  these  results. 

2.  Short  warning  times 

Ahrens  (1994,  p.  901),  Meissinger  (1995),  and  Solem  (1994,  p.  1015)  all  assume 
rectilinear  motion  for  short  warning  time  calculations.  Once  again,  this  approximation 
neglects  the  Earth's  gravitational  effects,  which  will  cause  the  NEO  to  deflect  towards  the 
Earth  in  a  hyperbolic  orbit. 

A  brief  look  at  hyperbolic  orbits  will  give  an  idea  of  how  Earth's  gravity  will  affect 
rectilinear  approximations.  When  a  sideways  motion  is  imparted  such  that  a 
(gravitationally  free)  miss  distance  of  1  R^  is  achieved,  the  result  is  to  establish  a 
hyperbolic  orbit  with  a  semiminor  axis  (b)  of  1  R^.  Using  a  typical  impact  velocity  of  20 
km/s,  actual  perigee  radius  (equation  from  Brown,  1992,  p.  27)  can  be  calculated  as 
follows: 


where 
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Voo  =  impact  velocity 
p  =  Gravitational  Constant  of  the  Earth 

Thus,  preliminary  approximations  of  energy  required  to  deflect  a  threat  may  result 
in  errors  of  14%  in  miss  distance,  in  the  terminal  case. 

This  thesis  intends  to  further  investigate  the  Earth's  gravitational  effects  on  energy 
required  to  defend  against  an  impacting  NEO. 
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m.  PROBLEM  FORMULATION 


A.  PROBLEM  STATEMENT 

Given  a  NEO  with  an  orbit  that  confirms  an  impending  collision  with  the  Earth, 
and  given  that  orbital  deflection  is  the  chosen  mitigation  goal,  it  is  necessary  to  know  the 
AV  (both  magnitude  and  direction)  required  to  deflect  the  asteroid  into  a  safe  orbit. 
Determination  of  this  required  AV  must  include  multi-body  gravitational  effects.  A 
detailed  imderstanding  of  third-body  perturbations  will  help  to  refine  previous  two-body 
approximations. 

The  goal  of  this  thesis  is  to  investigate  the  effects  of  the  Earth's  gravitation  on  the 
AV  required  to  deflect  an  Earth-impacting  NEO. 

B.  ASSUMPTIONS 

Several  assumptions  were  made  to  simplify  the  analysis. 

•  The  Earth  and  the  threat-NEO  are  in  co-planar  orbits. 

•  The  Earth  is  considered  to  be  in  a  perfectly  circular  heliocentric  orbit  at  a  radius 
of  1  AU. 

•  The  threat-NEO  is  originally  in  an  elliptical  heliocentric  orbit.  NEOs  in 
hyperbolic  or  parabolic  heliocentric  orbits  were  not  considered. 

•  The  control  maneuvers  are  impulsive. 

•  Other  than  the  impulse,  no  non-gravitational  forces  (such  as  outflowing  of  gas 
and  dust  in  comets)  are  included. 

C.  HYPERBOLIC  MAPPING  AND  ANALYSIS 

The  first  step  was  to  understand  how  a  heliocentric  elliptical  orbit  maps  into  a 
geocentric  hyperbolic  orbit.  The  primary  goal  was  to  determine  the  velocity  of  the  NEO 
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with  respect  to  the  Earth  (V^g^)  at  Earth  impact.  This  required  some  vector  analysis 

combined  with  orbital  mechanics  and  geometry. 

Once  an  equation  for  was  found,  an  understanding  of  the  dynamics  can  be 

aided  by  finding  its  maximum.  The  maximum  was  determined  in  two  ways.  First, 
Mathworks'  MATLAB®  was  used  to  numerically  create  a  3-D  plot  of  versus  the 

NEO's  heliocentric  semimajor  axis  (a )  and  its  heliocentric  eccentricity  (e ).  To  confirm 
the  numerical  result,  Warterloo's  symbolic  manipulator.  Maple  V™,  was  used.  The 
original  equation  for  was  maximized  analytically.  The  plot  was  identical  to  that  of 

the  MATLAB®  analysis,  thus  verilying  both  analyses. 

D.  PATCHED  CONIC  APPROXIMATION 

Once  the  hyperbolic  mapping  was  realized,  the  analysis  was  extended  for  use  in  a 
patched  conic  approximation.  The  advantage  of  using  this  approximation  was  that 
previous  two-body  analysis  could  be  adapted  to  accoimt  for  third  body  effects. 

The  most  useful  aspect  of  the  patched  conic  approximation,  for  this  analysis,  was 
the  B-plane  and  the  impact  radius  (^,).  Figure  3.1  illustrates  these  parameters.  The  B- 
plane  is  orthogonal  to  the  asymptote  of  the  approach  hyperbola  and  placed  at  a  large 
distance  fi'om  the  Earth.  The  impact  radius  (bj)  is  the  semiminor  axis  of  a  hyperbola  that 
has  a  perigee  equal  to  the  radius  of  the  surface  of  the  Earth. 
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Figure  3.1  B-plane  and  Impact  Radius  (from  Brown,  1992,  p.  117) 

The  impact  radius  is  a  function  of  ^nbo/^  » which  was  solved  for  in  the  hyperbolic 

analysis.  Once  this  radius  is  determined,  it  can  replace  previous  miss  distances  of  1  R®  in 
two-body  analysis.  Thus,  third-body  effects  are  accounted  for  while  maintaining  the 
relative  simplicity  of  two-body  analysis. 

EL  NUMERICAL  TRAJECTORY  OPTIMIZATION 

Park,  Elder,  and  Ross  (Park,  1997)  use  MATLAB®  's  sequential  quadratic 
programming  method  to  numerically  solve  a  constrained  optimization  problem  to 
determine  the  minimum  AV  for  deflecting  an  Earth-crossing  asteroid.  As  referred  to 
earlier,  this  code  used  a  two-body  approximation.  It  is  possible  to  include  the  effects  of 
third-body  gravitation  in  this  code  by  adding  the  impact  parameter  in  two  different  ways. 

Discussing  these  modifications  first  requires  a  brief  description  of  the  code.  The 
code  requires  a  target  separation  distance.  This  target  separation  is  the  minimum 
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allowable  separation  between  the  Earth  and  the  perturbed  NEO.  Following  an  impulse, 
the  Earth's  orbit  and  the  NEO's  orbit  are  propagated  using  the  two-body  equations  of 
motion  to  determine  the  minimum  separation.  The  optimization  problem  minimizes  the 
impulse  by  equalizing  the  propagated  and  the  target  separations. 

The  first  modification  of  the  code  is  to  increase  the  target  separation  distance.  If 
this  target  separation  is  increased  to  a  minimum  impact  radius  in  the  B-plane,  then  the 
Earth's  gravity  will  not  cause  an  impact. 

Second,  the  constraints  of  the  optimization  problem  can  be  modified.  In  its 
original  form,  the  constraints  require  that  the  minimum  separation  equal  some  fixed 
distance.  With  the  first  modification,  this  fixed  distance  is  the  impact  radius  of  the  NEO's 
original  orbit.  Once  the  NEO's  orbit  changes  (after  application  of  the  AV),  there  is  a  new 
impact  radius.  If  this  new  impact  radius  is  larger  than  the  original  impact  radius,  the  NEO 
will  still  impact  the  Earth.  The  solution  is  to  modify  the  constraint  to  include  the  impact 
radius  directly.  This  modification  requires  that  the  minimum  separation  be  equal  to  the 
impact  radius. 
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IV.  HYPERBOLIC  ORBIT  ANALYSIS 


A.  CONIC  SECTIONS 

Orbits  in  an  inverse  square  gravitational  field  take  the  shape  of  conic  sections: 
circular,  elliptical,  parabolic,  or  hyperbolic.  Circular  orbits  are  a  special  case  of  the 
elliptical  orbit.  These  two  types  of  orbits  are  defined  when  the  object's  energy  is 
insufficient  to  escape  the  gravitational  attraction  of  the  central  mass.  The  parabolic  orbit 
is  a  transition  between  elliptical  and  hyperbolic  orbits.  An  object  passing  a  central  mass 
on  a  parabolic  trajectory  would  reach  an  infinite  distance  from  the  mass,  but  with  zero 
velocity.  Any  object  with  a  velocity  less  than  that  of  a  parabolic  orbit  will  not  escape  the 
central  mass,  and  thus  will  be  in  an  elliptical  orbit.  Conversely,  an  object  with  a  greater 
energy  than  that  of  a  parabolic  orbit  will  escape  the  central  mass.  When  the  object  is 
traveling  faster  than  this  escape  velocity,  it  is  in  a  hyperbolic  orbit. 

B.  ELLIPTICAL  ORBITS 

Before  they  become  a  threat  to  Earth,  nearly  all  NEOs  are  established  in  an 
elliptical  orbit  about  the  Sun.  Figure  4.1  shows  the  geometry  of  an  elliptical  orbit. 


Figure  4.1  Elliptical  Orbit  (after  Brown,  1992,  p.  9) 

For  the  two-dimensional  case  considered  in  this  thesis,  only  the  semi-major  axis 
(a),  eccentricity  (e),  and  true  anomaly  (u)  are  required  to  define  an  orbit  and  position. 
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a  =  semimajor  axis,  the  distance  from  the  center  of  the  ellipse  to  the  long  edge 
e  =  eccentricity,  defines  the  shape  of  the  orbit  (0  <  e  <  1  for  an  elliptical  orbit) 

V  =  true  anomaly,  defines  the  position  within  the  orbit 

The  most  important  elliptical  orbit  relation  required  for  this  analysis  is  the  NEC's 
heliocentric  velocity  as  it  reaches  the  orbit  of  the  Earth.  This  velocity  is  given  by 
(Brown,  1992,  p.  19) 


Also  useful  is  the  distance  of  the  NEC  from  the 

anomaly,  given  by  (Brown,  1992,  p.  18) 

a[\-e^) 

r  - - 

1  +  ecosv 


(4.1) 

Sun  as  a  function  of  the  true 


(4.2) 


C.  HYPERBOLIC  ORBITS 

In  an  impact  scenario,  a  NEC  which  was  in  an  elliptical  orbit  about  the  Sun  will 
transition,  as  it  approaches  the  Earth,  to  a  hyperbolic  orbit  about  the  Earth.  An 
understanding  of  hyperbolic  orbits,  and  this  transition,  is  therefore  fundamental  to  studies 
of  the  orbital  mechanics  of  NEC  impacts. 

The  orbital  parameters  of  a  hyperbolic  orbit,  shown  graphically  in  Figure  4.2,  are 
similar  to  the  familiar  elliptical  orbit  parameters  (Brown,  1992,  p.  21): 

Tp  =  periapsis  radius 

a  =  absolute  magnitude  of  the  semimajor  axis,  the  distance  from  the  asymptote 
focus  to  the  periapsis 

b  =  semiminor  axis,  distance  from  the  asymptote  to  a  parallel  passing  through  the  central  body 
c  =  the  distance  from  the  asymptote  focus  to  the  center  of  mass 
e  =  eccentricity  =  ^  (greater  than  1  for  hyperbolic  orbit) 
jS  =  angle  of  the  asymptote 
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dg  =  trae  anomaly  of  the  asymptote 


Figure  4.2  Elements  of  a  hyperbola  (after  Brown,  1992,  p.  22) 


All  of  the  orbital  parameters  of  a  hyperbolic  orbit  can  be  determined  with  a 
knowledge  of  two  parameters:  the  velocity  at  an  infinite  distance  from  the  central  mass 
(Voo)  and  the  semiminor  axis  (b).  From  knowledge  of  V®,  the  semimajor  axis  can  be 
calculated  by  the  relation  (Brown,  1992,  p.  28) 


(4.3) 


Once  the  semimajor  axis  (a)  is  determined,  the  eccentricity  can  be  calculated  by 
(Brown,  1992,  p.  27) 


Periapsis  radius  (rp )  is  expressed  as  (Brown,  1992,  p.  27) 


(4.4) 
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And  finally,  the  angle  of  the  asymptote  (j8)  is  (Brown,  1992,  p.  27) 

tan  jS  =  — 
a 


D.  PATCHED  CONIC  METHOD 

The  patched  conic  method  was  originally  developed  to  simplify  the  planning  of 
interplanetary  trajectories,  yet  it  applies  directly  to  analysis  of  NEO  impacts  on  the 
Earth. 

Each  celestial  body,  due  to  its  mass,  generates  a  gravitational  field  which  affects  all 
other  celestial  bodies  according  to  Newton's  Laws  of  motion  and  gravity  (Wiesel,  1997,  p. 


=2.~3 - 


j*i 

This  is  known  as  the  N-body  problem.  An  exact  solution  for  more  than  two 
bodies  does  not  exist  (Weisel,  1997,  p.  35).  However,  the  patched  conic  method  makes 
some  minor  approximations  that  reduce  the  N-body  problem  to  solvable  two-body 
problems. 

The  patched  conic  method  assumes  that  an  object  is  influenced  by  the 
gravitational  field  of  a  planet  only  when  it  is  within  the  planet's  "sphere  of  influence." 
Beyond  the  sphere  of  influence,  the  object  is  considered  to  only  be  affected  by  the  Sun's 
gravitation.  The  radius  of  the  sphere  of  influence  is  somewhat  nebulous,  but  Laplace 
suggests  (Brown,  1992,  p.  97), 


^SOI  ~  ^ 
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where 


=  radius  of  the  sphere  of  influence  of  the  planet 
R  =  mean  orbital  radius  of  the  planet 

^planet  ~  planCt 

M^un  ~  mass  of  the  Sun 

For  Earth,  this  equates  to  a  sphere  of  influence  of  0.9  x  lO^km. 

In  the  case  of  a  NEO  impact  scenario,  the  NEO  begins  in  an  elliptical  orbit  about 
the  Sun.  Once  within  the  sphere  of  influence  of  the  Earth,  the  NEO's  motion  is  described 
by  two-body  orbit  equations  for  a  hyperbolic  orbit  about  the  Earth. 


E.  DETERMINING  V^o 


To  define  the  orbit  of  the  NEO  as  it  approaches  the  Earth,  it  is  necessary  to 
determine  Voo.  Figure  4.3  shows  the  vector  geometry  of  an  impact  scenario,  where: 

^NBo/  ~  velocity  of  the  NEO  with  respect  to  the  Sun 

Vffi/  =  the  velocity  of  the  Earth  with  respect  to  the  Sun 

/© 

~  velocity  of  the  NEO  with  respect  to  the  Earth,  defined  as  Voo  in  this 
thesis 

Strictly  speaking,  ^  V„ ,  because  the  NEO  is  not  at  infinity  with  respect  to 
the  Earth.  However,  in  the  patched  conic  method  the  approximation  is  made  as 


V  +  V  =  V 


% 


Using  vector  algebra,  Vo,  can  be  expressed  as: 

(4.9) 

V^=V.  =  V,„/^-Vy^  (4.10) 

The  velocity  of  the  Earth  with  respect  to  the  Sun  (V*,  )  can  be  determined  at  any 

7® 

time,  using  the  equation  for  velocity  in  a  circular  orbit 

v^=  [^X  (4.11) 


with 


W  “  distance  of  the  Earth  from  the  Sun  =  1  AU 

/® 


To  find  the  other  unknown. 


some  approximations  are  required.  The 


radius  of  the  sphere  of  influence  is  very  small  in  comparison  to  the  radius  of  the  Earth's 
orbit  about  the  Sun.  Thus,  calculating  the  velocity  of  the  NEO  as  it  crosses  the  Earth's 
orbit  is  a  very  close  approximation  to  the  velocity  as  it  crosses  the  Earth's  sphere  of 
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influence  (SOI).  In  the  worst  case  geometry,  this  approximation  introduces  an  error  in 
radius  of 


_  0.9  X 10  fern 


^xl00%  = 


/<!) 


1.496xl0«itm 


X 100%  =  0.6% 


(4.12) 


Strictly  speaking,  this  approximation  changes  equations  4.9  and  4.10  from 

equalities  to  approximations. 

Now  the  magnitude  of  calculated  by 

/® 


/  _  1^^® 

^NEO/  ~  ( 


(4.13) 


with 


a  =  the  semimajor  axis  of  the  NEO  with  respect  to  the  Sun 

For  the  direction  of  V^eq/  >  the  flight  path  angle  (y )  is  required.  The  flight  path 

/® 

angle  is  a  function  of  the  position  of  the  object  within  its  orbit.  This  position  is  defined 
by  the  true  anomaly  (u),  shown  in  Figure  4.4.  In  an  impact  scenario  with  a  known  NEO 
orbit,  the  true  anomaly  at  impact  can  be  calculated  as  follows  (where  e :  =  the  eccentricity 
of  the  NEO's  heliocentric  orbit): 

a(l-e^) 


NEOj 


0/ 

/® 


®  J  IMPACT 


=  V  = 

/© 


1+ecos  V 


'NEO. 


% 


(4.14) 


IMPACT 


After  manipulation.  Equation  4.14  yields 


^NEO/ 

1  =  cos'‘ 

\  7© 

J  IMPACT 

era 


©/ 

/® 


(4.15) 
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Figure  4.4  Definition  of  true  anomaly 
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Equations  4.18  are  an  important  result.  Although  they  look  complicated,  it  is 
apparent  that  ¥„  is  a  function  of  only  two  variables:  the  NEO's  heliocentric  eccentricity 
(e)  and  semimajor  axis  (a ).  The  flight  path  angle  and  radius  are  set  by  the  geometry  of 
the  impact.  Since  the  Earth  is  assumed  to  be  in  a  circular  orbit,  and  the  problem  is  planar, 
this  geometry  is  defined  by  the  NEO's  orbit. 

F.  MAXIMUMVoo 

As  previously  discussed,  ¥„  is  a  critical  value  in  defining  the  hyperbolic  orbit  of 
the  NEO  about  the  Earth's  center  of  gravity.  Now  that  V*  has  been  directly  related  to  e 
and  a ,  a  better  understanding  of  this  relation  can  be  achieved  by  attempting  to  find  its 
maximum. 


1. 


Variable  Bounds 


If  it  is  possible  to  find  bounds  for  each  parameter  in  equations  4. 16,  it  may  be 
possible  to  detamine  a  maximum  Vo,.  Such  a  study  will  potentially  provide  an  insight 
into  the  full  geometry  of  the  problem. 

To  determine  a  value  for  the  minimum  semimajor  axis  of  the  NEO  (amm),  recall  that 


a  = 


a  p 


(4.19) 


A  minimum  value  for  the  NEO's  semimajor  axis  is  mandated  by  the  fact  that  the 
orbit  must  intersect  the  Earth's  orbit  for  an  impact  to  occur.  This  sets  the  value  of  ra  to 
the  radius  of  the  Earth's  orbit.  The  NEO  also  cannot  have  a  perihelion  radius  less  than  the 
radius  of  the  Sxm  (Rg^.  Thus,  amin  is  defined  by 


^min 


(4.20) 
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hand,  impact  occurs  at  the  NEC's  aphelion,  =  tc  radians.  In  between  these 

V  J  IMPACT 

two  extremes,  the  impact  can  happen  at  one  of  two  separate  true  anomalies.  Figure  4.4 
(above)  demonstrates  this  fact.  Thus,  the  true  anomaly  at  impact  in  bounded  as  follows 
(expressed  in  radians): 

(4.24) 

To  understand  the  bounds  of  the  flight  path  angle  (y ),  a  graph  of  the  relation 
between  flight  path  angle  and  true  anomaly  is  useful.  This  graph  is  presented  in  Figure 
4.5. 
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Figure  4.5  Flight  Path  Angle  of  the  NEO  at  impact  relative  to  True  Anomaly  for  various 

eccentricities 


7t  71 

One  can  see  from  the  figure  that  the  flight  path  angle  varies  between - and  — 

Jit  Jt. 


71  71 

I - <  y  <  — 

'2^2 


(4.25) 


In  conclusion,  a  study  of  the  bounds  of  each  parameter  does  not  readily  reveal  a 
relation  for  ¥«>.  Thus,  a  plot  of  V*  may  be  more  helpful. 


2.  Graphical  Display 

Since  Voo  is  a  function  of  only  two  variables,  it  is  possible  to  generate  a  three-dimensional 
plot  of  Voo  against  e  and  a .  A  small  MATLAB  m-file  was  written  to  generate  such  a 
plot.  The  m-file  is  listed  in  Appendix  A.  The  plot  is  displayed  in  Figure  4.6. 
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VJnfinily  (km/s) . 


Figure  4.6  as  a  function  of  NEO  eccentricity  and  semimajor  axis 

The  flat  portion  of  this  graph  represents  NEO  orbits  that  do  not  intersect  the 
Earth's  orbit  Note  that  for  a  given  NEO  eccentricity  there  is  a  semimajor  axis  that 
generates  the  maximum  Vo,.  This  is  demonstrated  graphically  by  the  plot  in  Figure  4.7. 
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A  second  effect  apparent  in  Figure  4.7  is  the  increase  in  Vo,  as  the  NEO's 
eccentricity  increases. 

3.  Closed  Form  Solution 

As  a  final  attempt  to  understand  the  geometry  that  defines  V®,  the  symbolic 
manipulator  Maple  V™  was  used  to  generate  a  closed  form  solution  for  the  maximum 
value  for  Voo.  The  Maple  V™  code  is  included  in  Appendix  B. 

The  first  step  in  the  effort  to  find  a  closed  form  solution  is  to  combine  all  of 
Equations  4.18  into  one  equation.  Because  of  the  assumption  that  the  NEO  is  in  an 
elliptical  orbit  about  the  Sun,  it  is  possible  to  reduce  the  Maple  V™  solution  to  the 
following  equation: 
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(4.26) 

Next,  the  square  of  the  magnitude  of  the  velocity  was  calculated  by  adding  the 
square  of  both  components  of  the  vector  in  Equation  4.26.  In  order  to  simplify  the 
expression  for  the  location  of  maximum  Va,,  the  function  was  left  as  the  square  of  the 
magnitude.  By  eliminating  the  square  root,  the  expression  is  simplified,  but  the  final 
result  will  be  the  same.  After  simplification,  the  magnitude  of  V®  reduces  to 


(4.27) 

Since  the  plots  in  the  last  section  show  that  there  is  a  semunajor  axis  (a)  for  each 
eccentricity  (e)  which  generates  a  maximum  Voo,  it  will  be  most  instructive  to  find  an 
expression  for  maximum  Vo,  with  a  fixed  eccentricity.  Thus,  Equation  4.27  is 
differentiated  with  respect  to  the  semimajor  axis,  with  the  eccentricity  held  constant. 
This  derivative  reduces  to 
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(4.28) 

Solving  Equation  4.28  for  the  semimajor  axis  gives 

(4.29) 

This  relatively  simple  equation,  considering  the  complexity  of  the  derivation,  is  an 
important  result.  The  equation  determines  the  semimajor  axis  that  will  generate  the 
maximum  Voo  for  a  given  eccentricity. 

Mathematically,  there  is  no  guarantee  that  Equation  4.29  is  a  maximum.  Since 
only  one  derivative  has  been  performed,  all  that  is  known  is  that  Equation  4.29  is  an 
extremum.  In  the  Maple  V™  procedure,  inputting  Equation  4.29  into  the  second 
derivative  produced  a  very  complicated  result.  Analysis  of  this  result  did  not  confirm 
that  Equation  4.29  was  a  maximum.  However,  an  eccentricity  of  0.9  gave  a  semimajor 
axis  of  1 .739,  which  matches  the  maximum  of  the  plot  in  Figure  4.7.  This  confirms  that 
Equation  4.29  matches  the  graphical  result  and  yields  a  maximum.. 


G.  MAPPING  TO  A  GEOCENTRIC  ORBIT 

As  stated  before,  two  parameters  are  required  to  define  a  hyperbolic  orbit:  Voo  and 
the  semiminor  axis  (b).  The  last  sections  went  into  great  detail  about  determining  Voo- 
However,  knowledge  of  the  NEO's  heliocentric  orbit  is  insufficient  to  uniquely  fix  the 
geocentric  semiminor  axis. 

In  order  to  fix  the  semiminor  axis,  one  would  require  exact  knowledge  of  the 
location  of  the  NEO,  with  respect  to  the  Earth,  as  it  crosses  the  Earth's  sphere  of 
influence.  All  that  is  known  about  the  impact  problem  is  the  NEO's  heliocentric  orbital 
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parameters,  and  the  prediction  of  a  future  impact  with  Earth.  Knowledge  of  the 
heliocentric  orbit  only  fixes  the  NEO's  Voo  as  it  crosses  the  sphere  of  influence.  However, 
there  is  no  single  semiminor  axis  that  will  define  an  impact  based  on  ¥«,. 

A  parameter  used  in  the  patched  conic  approximation  for  interplanetary  arrival 
targeting  can  be  applied  directly  to  the  impact  problem.  This  parameter  will  determine  a 
minimum  semrmajor  axis  required  to  avoid  impact.  The  definition  of  fiiis  parameter  is 
first  based  on  the  B-plane.  Shown  in  Figure  3.1,  the  B-plane  is  perpendicular  to  the 
asymptote  of  the  approach  hyperbola,  and  placed  at  an  infinite  distance  from  the  Earth 
(Brown,  1992,  p.  1 16). 

In  this  Figure  the  semiminor  axis  (b)  is  a  radius,  drawn  in  the  B-plane,  from  the 
location  that  the  NEO  pierces  the  B-plane  to  a  line  perpendicular  to  the  B-plane  that 
intersects  the  center  of  mass  of  the  Earth.  The  impact  radius  (bi)  is  the  semiminor  axis 
that  defines  a  hyperbolic  orbit  which  is  tangent  to  the  surface  of  the  Earth.  Whenever  the 
NEO  passes  inside  the  impact  radius,  it  will  impact  the  Earth.  Thus,  the  impact  radius 
defines  a  minimum  safe  semiminor  axis.  The  impact  radius  is  a  function  of  Voo  and  can  be 
calculated  from  the  general  hyperbolic  relation  (Brown,  1992,  p.  28) 

For  the  impact  radius  specifically  about  the  Earth, 

=  1-^  +  1  (4-31) 

Since  the  previous  analysis  determined  Voo  as  a  function  of  the  NEO's  heliocentric 
orbital  parameters,  it  is  now  possible  to  define  the  Impact  Radius  as  a  function  of  the 
NEO's  orbital  parameters.  A  plot  of  this  relation  is  included  in  Figure  4.8. 
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Figure  4.8  Impact  Radius  as  a  function  of  the  NEO's  eccentricity  and  semimajor  axis 


Note  that,  for  clarity,  the  axes  have  been  rotated  from  those  used  in  Figure  4.6. 

Figure  4.8  confirms  an  important  conclusion  that  can  be  reached  from  study  of 
Equation  4.31:  the  impact  radius  increases  significantly  as  ¥«>  decreases.  This  is  an 
important  result  when  extended  to  analysis  of  the  NEO  impact  problem.  Hence,  the 
greatest  error  in  previous  analysis  (Knudson,  1995;  Park,  1997;  Elder,  1997;  Ahrens 
1994;  Meissinger,  1995;  Solem,  1994)  which  neglected  third-body  effects  arises  when 
the  NEO's  orbit  defines  a  small  V<x>.  These  orbits  are  nearly  circular,  and  close  to  the 
Earth  (that  is,  with  a  ~  1  AU). 
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V. 


APPLICATION  OF  THE  PATCHED  CONIC  APPROXIMATION 


The  preceding  chapter  established  a  theoretical  foundation  for  applying  a  patched 
conic  approximation  to  a  three-body  analysis  of  NEO  hazard  mitigation.  Using  this 
approach,  a  two-body  analysis  can  be  modified  to  include  third-body  effects. 

In  their  paper  entitled  "Minimum  Delta-V  For  Deflecting  Earth-Crossing 
Asteroids,"  Park,  Elder,  and  Ross  (Park,  1997)  presented  a  trajectory  optimization 
problem  that  was  solved  numerically  using  a  sequential  quadratic  programming  (SQP) 
method.  The  goal  of  the  code  was  to  determine  an  optimal  AV  to  deflect  an  asteroid,  but 
it  did  so  using  only  a  two-body  approximation.  The  patched  conic  approximation  has 
since  been  added  to  the  code.  The  results  are  presented  here. 


A.  ORIGINAL  TRAJECTORY  OPTIMIZATION  PROBLEM 


The  performance  index  chosen  to  minimize  AV  is 
J=^AV^^  +  AV^ 


(5.1) 


where  AV||  is  the  component  of  AV  which  is  parallel  to  the  motion  of  the  object  and  AVj^ 
is  the  component  perpendicular  to  the  motion.  The  constraints  consist  of  the  two-body 
equations  that  govern  the  motion  of  the  NEO  and  the  Earth  about  the  Sim.  Three  terminal 
boundary  conditions  are  set  by  the  desired  minimum  miss  distance.  First,  the  miss 
distance  at  the  final  time  must  be  equal  to  the  minimum  miss  distance.  Expressed 
mathematically 


(5.2) 


where  R  is  the  distance  between  the  Earth  and  the  NEO  and  Rcriticai  is  the  proposed  miss 
distance  firom  the  Earth.  The  two  additional  boundary  conditions  come  from  the 
requirement  that  R  be  at  a  minimum  when  it  reaches  Rcriticai-  Because  R  is  continous  and 
differentiable,  this  can  be  expressed  as 
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R  =  0 


(5.3) 

(5.4) 


R>0 

B.  TfflRD  BODY  MODIFICATIONS  TO  THE  ORIGINAL  PROBLEM 

Third  body  effects  can  be  added  to  the  original  formulation  by  modifying  the 
boundary  conditions  in  one  of  two  ways.  Both  modifications  are  new  definitions  for 
Rcriticai  as  used  in  Equation  5.2. 

In  the  first  modification,  Rcnticai  is  redefined  to  be  equal  to  the  impact  radius  (in 
the  B-plane)  of  the  original  NEO  orbit.  Once  the  original  orbital  parameters  of  the 
impacting  NEO  are  determined,  the  user  can  input  the  desired  miss  distance  and  the 
modified  code  computes  the  impact  radius  as  a  function  of  the  original  NEO  parameters. 
This  can  be  expressed  mathematically  as 

R  —  separation^  =  0  (5-5) 

where  "separation"  is  the  desired  minimum  separation  between  the  NEO  and  the  center  of 
mass  of  the  Earth. 

The  second  modification  defines  Rcnticai  as  the  impact  radius  of  the  NEO  orbit 
after  perturbation.  Instead  of  defining  Rcnticai  based  on  the  original  NEO  orbital 
parameters,  this  modification  recalculates  the  impact  radius  after  the  perturbation  is 
applied  and  confirms  that  Rcnticai  is  still  equal  to  the  impact  radius.  The  second 
modification  therefore  changes  the  boundary  condition  to 

R  -  bj{a,e,  separation)  =  0  (5.6) 

The  two  modifications  produced  slightly  different  results.  The  first  modification 
changed  only  one  calculation  of  bj  and  thus  had  little  effect  on  the  run  time  for  the  code. 
The  second  modification  required  several  calculations  of  bj  during  execution  of  the  code 
and  thus  increased  the  run  time  significantly.  However,  Equation  5.6  is  a  more  direct 
statement  of  the  required  miss  distance.  This  equation  accounts  for  the  fact  that  the 
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impulse  which  diverts  the  NEO  will  generate  a  new  orbit.  Since  the  impact  radius  is  a 
function  of  the  NEO's  orbit,  the  perturbation  will  also  generate  a  new  impact  radius. 
Equation  5.6  thus  ensures  that  even  the  perturbed  NEO  will  not  pass  within  the  impact 
radius.  All  results  discussed  below  refer  to  the  third-body  modification  expressed  in 
Equation  5.6. 


G.  SQP  NON-LINEAR  OPTIMIZATION  ALGORITHM 


The  trajectory  optimization  code  written  by  Dr.  Park  is  called  Earth  Defense 
against  Asteroid  Impact  (EDAJ). 

From  the  user's  viewpoint,  the  code  flows  according  to  the  diagram  in  Figure  5.1. 


Figure  5.1  EDAI  Flow  Diagram 
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At  this  early  stage  in  the  development  of  EDAI,  the  user  is  required  to  be  veiy 
knowledgeable  about  the  code's  design.  There  are  several  parameters  that  must  be 
adjusted,  within  the  code,  before  each  run.  With  some  adjustments,  discussed  in  the 
concluding  chapter,  the  code  can  be  much  more  user-friendly. 

As  with  all  numerical  optimization  programs,  the  EDAI  code  is  somewhat 
unstable  with  bad  initial  guesses.  A  critical  parameter  in  the  code  is  the  initial  guess  of  the 
minimum  AV.  If  this  initial  guess  is  significantly  different  from  the  true  value,  the 
numerical  optimization  can  become  unstable^  as  is  typical  with  numerical  codes. 

Currently,  there  is  no  way  to  change  this  initial  guess  without  getting  into  the  code. 

Thus,  if  the  user  wishes  to  initialize  a  new  EDAI  run  with  NEO  orbital  parameters  that 
are  different  from  the  previous  run,  he  must  input  an  accurate  initial  guess  into  the  code. 
The  initial  guesses  used  for  this  analysis  are  included  in  Appendix  C. 

Another  instability  comes  from  round-off  error  when  trying  to  achieve  a  minimum 
separation  of  one  Earth  radius.  Since  the  code  uses  AU  as  the  distance  unit,  a  minimxun 
separation  of  1  is  extremely  small.  This  is  easily  solved,  due  to  the  linearity  of  the 
problem,  by  running  the  code  vsdth  a  larger  minimum  separation  (e.g.,  1 0  and  dividing 

the  result  to  achieve  the  desired  separation  (e.g.,  divide  by  10). 

The  results  of  each  run  are  saved  in  a  file  which  is  named  within  the  EDAI  code. 
This  technique  allows  for  many  kinds  of  analysis  of  the  results.  However,  in  order  to 
change  the  file  name  one  must  make  a  modification  directly  to  the  code. 

D.  RESULTS  OF  THIRD  BODY  MODIFICATIONS 

The  impact  radius  (bj)  is  highly  dependent  on  the  orbital  parameters  of  the  NEO, 
as  displayed  in  Figure  4.9.  Since  the  third  body  modifications  to  the  EDAI  code  make  use 
of  the  impact  radius,  the  differences  in  results  are  also  highly  dependent  on  the  NEO's 
orbital  parameters. 
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First,  consider  NEOs  with  a  nearly  circular  orbit  (e  close  to  zero).  In  this  case,  the 
impact  radius  is  relatively  large.  Thus,  one  would  expect  the  difference  between  the 
results  of  the  2-body  model  and  the  3 -body  model  to  be  relatively  large.  Figure  5.2  is  a 
plot  of  the  minimum  AV  required  to  generate  a  1  miss,  as  calculated  by  both  the  2- 
body  and  3 -body  codes,  for  a  nearly  circular  orbit.  (Note  that,  for  standardization,  all 
results  will  be  presented  at  1  minimum  separation.)  The  x-axis  is  the  time  that  the 
impulse  is  applied,  in  units  of  NEO  periods  before  impact  with  the  Earth.  The  y-axis  is 
the  minimum  AV  required  in  units  of  cm/s. 

Park  (Park,  1997,  p.  6)  presents  a  thorough  explanation  of  the  results  of  the  EDAI 
code.  In  summary,  there  are  two  effects  apparent  in  the  results.  The  first  order  variation 
in  minimum  AV  required  is  dependent  on  the  change  in  orbital  elements  of  the  NEO  after 
perturbation.  The  second  order  variation  shows  local  minima  at  the  perihelion  of  the 
NEO  orbit.  These  extrema  demonstrate  the  commonly  known  effect  that  a  AV  applied  at 
the  perihelion  has  a  greater  efficiency  in  changing  an  orbit  than  when  applied  at  any  other 
place. 
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Figure  5.2  Minimum  AV  for  1  Separation  (e=0.1,  a=1.05  AU) 

For  the  original  NEO  orbit  used  in  the  plot  for  Figure  5.2,  the  impact  radius  is 
4. 192  R^.  Figure  5.3  shows  that  when  the  results  of  the  three-body  analysis  are  divided 
by  the  results  of  the  two-body  analysis  the  answer  oscillates  about  a  factor  of 
approximately  4. 195.  This  oscillation  could  be  identical  to  the  variation  in  the  impact 
radius  for  various  impulse  times.  However,  this  hypothesis  has  not  yet  been  verified. 
Verification  is  suggested  in  the  recommendations  for  future  work. 


44 


i  Y(3-bodyV/kY(2-body)  vs.  Impulse  Time 


Figure  5.3  AV(3-body)/AV(2-body)  (e=0.1,  a=1.05  AU) 

Figure  5.4  plots  the  difference  between  the  3-body  and  the  2-body  models.  The 
critical  information  in  this  plot  is  the  increasing  effect  of  the  Earth's  gravity  as  the  impulse 
time  decreases.  This  demonstrates  the  need  to  include  the  Earth's  gravity  in  planning  for 
short  warning  deflection  scenarios.  By  not  including  this  effect,  the  required  minimum  AV 
can  be  in  error  by  more  than  an  order  of  magnitude. 


45 


Figure  5.4  Difference  Between  3-Body  and  2-Body  models  (e=0.1,  a=l  .05  AU) 


Referring  once  again  to  Figure  4.9,  it  is  apparent  that  the  impact  radius  for  highly 
eccentric  NEOs  is  nearly  one.  Thus,  one  would  expect  that  adding  third  body  effects  to 
the  analysis  of  such  NEOs  would  have  little  effect.  Figure  5.5  is  a  plot  of  the  minimum 
AV  required  for  a  NEO  with  eccentricity  increased  to  e=0.9  and  the  semi-major  axis 
remaining  at  a=l  .05  AU.  While  it  is  nearly  impossible  to  tell,  both  the  2-body  and  3- 
body  results  are  included  in  the  plot.  This  very  close  correlation  is  a  result  of  the  fact 
that  the  impact  radius  for  this  original  NEO  orbit  is  1 .0593  R^. 
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Figure  5.5  Minimum  AV  for  1  Separation  (e=0.9,  a=l  .05  AU) 

While  the  effects  of  the  third-body  modification  are  significantly  smaller  than  with 
the  nearly  circular  orbit,  the  trends  are  the  same  with  both  NEO  orbits.  Figure  5.6  shows 
that  the  error  increases  with  decreasing  impiilse  time,  just  as  in  the  first  orbit. 
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impulse  Time  (periods] 


Figure  5.6  Difference  Between  3-Body  and  2-Body  models  (e=0.9,  a=1.05  AU) 

A  critical  benefit  from  adding  the  third  body  effects  to  the  code  is  that  the  short 
warning  results  are  now  much  more  accurate.  The  final  two  plots  are  details  of  the 
minimum  AV  results  for  the  two  orbits  for  the  timeframe  of  less  than  one  orbital  period. 
Note  the  sharp  variation  in  the  plot  in  Figure  5.8.  The  variation  occurs  when  the  NEO  is 
near  perihelion.  As  in  the  subsequent  passes  through  perihelion,  the  AV  requirement  to 
perturb  the  NEO  into  a  safe  orbit  is  minimum  at  perihelion.  Figure  5.8  shows  that  this 
applies  during  short  warning  time  scenarios  as  well. 
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Figure  5.7  Short  Warning  Minimum  AV  for  1  Separation  (e=0.1,  a=1.05  AU) 
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VL  MITIGATION  CAPABILITY  OF  CURRENT  IMPULSE  TECHNOLOGIES 


With  a  high  degree  of  confidence  in  the  new  three-body  minimum  AV  results,  it  is 
now  possible  to  evaluate  the  potential  of  current  technologies  for  deflection  of  NEOs. 
While  Chapter  II  presented  several  ideas  with  potential  for  NEO  deflection,  the  three- 
body  analysis  of  this  thesis  was  based  on  an  assumption  of  an  instantaneous  AV  (that  is, 
an  impulse).  An  impulse  is  practically  impossible,  and  the  only  three  methods  discussed 
which  approximate  an  impulse  are  kinetic  impactors,  a  stand-off  nuclear  blast,  and  a 
surface  nuclear  blast. 

There  are  many  variables  to  consider  when  analyzing  the  capability  of  an  impulse 
technology,  including  the  physical  properties  of  the  NEO,  the  mass  and  closing  velocity 
of  a  kinetic  impactor,  or  the  explosive  yield  and  stand-off  distance  of  a  nuclear 
interceptor.  Ahrens  and  Harris  (Ahrens,  1994,  p.  897)  present  a  detailed  analysis  of  the 
effects  of  these  variables.  This  analysis  was  combined  with  the  minimum  AV  results  from 
the  last  chapter  to  evaluate  current  deflection  capabilities. 

A.  KINETIC  IMPACT 


A  kinetic  impactor  imparts  a  AV  on  a  NEO  by  the  exchange  of  momentum  from 
the  collision  of  the  impactor's  mass  and  by  creating  a  crater.  The  ejecta  from  the  crater 
accelerates  away  from  the  NEO,  thus  causing  a  change  in  momentum. 

In  their  analysis,  Ahrens  and  Harris  obtain  an  expression  for  the  mass  ratio 
between  the  kinetic  impactor  (Mj)  and  the  NEO  (Mneo)  as  follows.  (Ahrens,  1994,  p. 


906) 


M, 


M, 


mo 


1+0.16 


VPi) 


0309 


(6.1) 


where 

Av  =  change  in  NEO  velocity 
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V.  =  impact  speed 
p  =  NEO  density 
p.  =  kinetic  impactor  density 
Y  =  NEO  material  strength 

The  numerator  of  Equation  6.1  gives  the  AV  resulting  from  an  inelastic  collision, 
while  the  denominator  accounts  for  the  added  momentum  change  from  the  crater  ejecta. 

It  should  be  noted  that  Ahrens  and  Harris  differentiate  between  two  different 
regimes  in  their  analysis.  For  smaller  NEOs,  the  strength  of  the  NEO  dominates  the 
cratering  mechanism.  This  effect  is  referred  to  as  the  strength  regime.  The  gravity  regime 
applies  to  larger  NEOs  where  gravitational  effects  dominate  cratering  (Ahrens,  1994,  p. 
904).  The  results  were  not  significantly  different,  so  the  strength  regime  is  used  for  the 
current  analysis. 


B.  STAND-OFF  NUCLEAR  BLAST 


A  stand-off  nuclear  blast  can  impart  a  AV  on  a  NEO  by  heating  its  surface  imtil  it 
expands  and  ablates  away  from  the  NEO.  The  ablated  material  is  mass  accelerating  away 
from  the  NEO,  thus  causing  a  momentum  change.  Ahrens  and  Harris  (Ahrens,  1994,  p. 
912)  conclude  that  the  effect  of  a  stand-off  nuclear  blast  can  be  expressed  as 


Av»0.1 


nAW 

If 


(6.2) 


with 

n  =  efficiency  of  neutron  production 

A  =  Geometrical  efficiency  factor  (accounts  for  stand-off  distance) 

W  =  nuclear  explosive  yield  (expressed  in  Kt) 

D  =  NEO  diameter  (in  km) 

Equation  6.1  expresses  the  NEO  size  in  terms  of  mass,  while  Equation  6.2  uses 
diameter.  For  consistency,  and  for  more  direct  application  in  the  MATLAB  m-file. 
Equation  6.2  was  modified  to  give 


52 


^NEO  ~ 


0.  \nAWp7t 
6Av 


(6.3) 


C.  SURFACE  NUCLEAR  BLAST 

The  mechanism  for  momentum  change  with  a  surface  nuclear  blast  is  nearly 
identical  to  that  of  the  kinetic  impactor,  except  that  the  increased  energy  of  the  nuclear 
explosive  creates  a  significantly  larger  crater.  Since  there  is  little  experimental  data  about 
the  effects  of  nuclear  cratering  in  the  strength  regime,  Ahrens  and  Harris  present  only  a 
rough  but  educated  estimate.  The  effect  of  a  surface  nuclear  blast  is  therefore 
approximated  by  (Ahrens,  1994,  p.  913) 

4e-5AvM^£o  (6.4) 

where  Av  is  expressed  in  cm/s  and  Mneo  in  kg- 

D.  ASSUMPTIONS 

Several  assumptions  are  required  in  the  analysis  of  mitigation  technology 
capabilities.  The  physical  properties  of  the  NEO  required  by  Equations  6.1  and  6.3  are 
the  yield  strength  (Y)  and  the  density.  The  yield  strength  can  vary  from  lO’  dyne/cm^  for 
soft  rock  or  ice  to  lO’  dyne/cm^  for  hard  rock  (Ahrens,  1994,  p.  906).  For  this  analysis, 
a  value  of  10*  dyne/cm^  was  used.  The  density  of  asteroids  ranges  from  2x10^  kg/m^  to 
5x10*  kg/m*  with  a  mean  density  of  3x10*  kg/m*  (Elder,  1997,  p.  10).  Comets  are 
somewhat  less  dense  than  asteroids,  with  estimates  of  1000  to  2000  kg/m*  (Elder,  1997, 
p.  10).  A  mean  asteroid  density  of  3x10*  l^m*  was  used  for  this  analysis. 

Assumptions  are  also  required  for  the  mitigation  technology.  The  goal  of  this 
analysis  is  to  determine  the  maximum  NEO  size  that  can  be  diverted  using  a  single 
mission  with  today's  technology.  It  is  always  possible  to  send  multiple  diversion 
missions  against  a  NEO  threat,  but  a  single  mission  baseline  gives  the  most  fundamental 
result.  Rustan  (Rustan,  1994,  p.  1070)  presents  a  summary  of  laimch  vehicle 
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performance  which  includes  the  useful  payload  that  can  be  place  on  an  interplanetary 
trajectory.  He  indicates  that  the  Russian  Energia  vehicle  can  place  a  payload  mass  of 
18,040  kg  into  an  inteiplanetary  (or  NEO  intercepting)  trajectory.  Equation  6.1  implies 
that  an  impactor  with  high  density  would  be  most  effective  in  diverting  a  NEO.  Lead  has 
high  density  and  is  inexpensive  enough  to  be  a  reasonable  candidate  for  use  as  a  kinetic 
impactor.  While  the  impactor  would  involve  both  guidance  systems  and  concentrated 
mass,  a  fair  approximation  is  to  use  the  density  of  lead,  which  is  1 1  g/cm^.  The  closing 
velocity  of  the  kinetic  impactor  is  a  function  of  the  trajectory  chosen  to  achieve  intercept. 
For  this  analysis,  a  median  closing  velocity  of  20  km/s  was  chosen. 

A  search  of  imclassified  documentation  determined  that  the  largest  single  nuclear 
warhead  yield  in  a  currently  deployed  strategic  system  is  24  Mt  (Janes,  1996).  This  is 
the  warhead  used  in  the  Russian  SS-18  Mod  1  missile,  appropriately  nicknamed  "Satan." 
Since  no  information  about  the  warhead  mass  was  available,  it  was  assumed  to  be  less 
than  the  1 8,040  kg  useful  payload  of  the  Energia  launch  vehicle.  The  neutron  production 
efficiency  (n)  can  vary  from  0.03  to  0.3  (Elder,  1997,  p.  43).  A  median  value  of  n  =  0.15 
was  chosen  for  this  analysis.  The  geometrical  efficiency  factor  (A)  is  taken  from  a 
standoff  distance  of  0.4  NEO  radii  to  give  a  value  of  A  =  0.3  (Ahrens,  1994,  p.  912). 

Finally,  the  results  assume  that  the  AY  is  applied  in  an  optimal  direction. 

All  assumptions  used  for  analysis  of  the  capability  of  current  mitigation 
technologies  are  summarized  as  follows: 
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NEO  Physical  Properties 

Yield  strength 

Y  =  10*  dyne/cm^ 

Density 

p  =  3x10^  kg/m^ 

Kinetic  Impactor 

Mass 

Mi  =18,040  kg 

Density 

Pi  =11  g/cm^ 

Closing  Velocity 

Vi  =  20  km/s 

Nuclear 

Yield 

W  =  24  Mt 

Neutron  Production 

n  =  0.15 

Geometrical  Efficiency 

A  =  0.3 

Optimal  Impulse  Direction 

Single  Diversion  Mission. 

Table  6.1  Assumptions  Used  in  Mitigation  Capability  Analysis 


E.  TOUTATIS 

The  asteroid  Toutatis  is  an  Earth-Crossing  Asteroid  that  passes  within  10  lunar 
distances  every  four  years  (Park,  1997,  p.  9).  Its  low  inclination  makes  it  suitable  for 
analysis  with  the  planar  methods  used  in  this  thesis.  The  orbital  parameters  of  Toutatis 
are  as  follows: 

a  =  2.5154  AU 
e  =  0.6361 
i  =  0.47° 

To  apply  the  analysis  to  Toutatis,  it  must  be  assumed  that  i=0°.  This  condition 
is  necessary  to  create  an  impact  with  the  Earth,  as  well  as  to  conform  with  the 
assumptions  made  in  the  patched  conic  analysis. 
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Toutatis  takes  the  shape  of  two  attached  spheres  which  combine  to  give  a 
diameter  of  approximately  4.3  km  (Elder,  1997,  p.  45). 

Figure  6.1  indicates  the  capability  of  the  three  impulse  technologies  against  the 
Toutatis-type  asteroid.  The  Figure  shows  a  plot  of  the  NEO's  diameter  versus  the 
impulse  time  in  Earth  years.  It  should  be  noted  that  the  critical  properties  that  define  the 
amount  of  energy  required  for  mitigation  is  the  minimum  AV  and  the  NEO  mass.  The 
NEO's  density  was  assumed  above.  This  assumption  allows  calculation  of  NEO 
diameter,  assuming  a  spherical  NEO.  Since  the  NEO  diameter  is  somewhat  easier  to 
identify  v^th  than  its  mass,  the  diameter  is  used  in  these  plots. 


Figure  6.1  Mitigation  Capability  Against  NEOs  in  Toutatis-type  Orbit 


It  is  apparent  fi'om  Figure  6.1  that  Toutatis-type  NEOs  caimot  be  deflected  with  a 
single  kinetic  impulse  for  many  years.  The  diameter  of  Toutatis  (4.3  km)  lies 
comfortably  within  the  capability  of  either  nuclear  option.  However,  Figure  6.2  shows 
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that  a  short  warning  impact  from  Toutatis  caimot  be  deflected  with  any  single  mission. 
This  Figure  shows  the  deflection  capability  within  one  period  of  Toutatis.  From  this 
Figure,  it  is  apparent  that  a  surface  nuclear  impulse  can  deflect  a  Toutatis-type  asteroid 
when  the  impulse  is  applied  as  close  as  2  Earth  years  before  impact.  If  the  impulse  is 


applied  any  closer  than  that,  multiple  missions  will  be  required. 


Figure  6.2  Short  Warning  Mitigation  Capability  Against  NEOs  in  a  Toutatis-like  Orbit 


F.  NEREUS 

For  comparison,  a  second  asteroid  in  a  more  circular  orbit  that  Toutatis  was 
analyzed.  The  orbital  parameters  of  the  asteroid  Nereus  meet  the  criteria,  and  are  listed 
below  (NASA  Ames): 

a  =  1.4897  AU 
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e  =  0.3602 


i=  1.41° 

D  =  0.8km 

Although  this  asteroid  is  significantly  smaller  than  Toutatis,  the  long  range  plot  in 


Figure  6.3  indicates  that  the  kinetic  impactor  is  still  insufficient  as  early  as  18  years 
before  impact.  However,  it  can  be  deflected  with  a  single  mission  as  close  as  0.2  Earth 
years  before  impact,  as  indicated  in  Figure  6.4. 


Figure  6.3  Mitigation  Capability  Against  NEOs  in  Nereus-type  Orbits 
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Figure  6.4  Short  Warning  Mitigation  Capability  Against  NEOs  in  a  Nereus-like  Orbit 
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m  SUMMARY  AND  RECOMMENDATIONS  FOR  FURTHER  WORK 


A.  SUMMARY 

The  gravitational  effects  of  the  Earth  have  the  strongest  influence  on  minimum  AV 
calculations  for  deflecting  NEOs  in  nearly  circular  heliocentric  orbits.  Incorporating  the 
Earth's  gravitational  effects  can  increase  the  minimum  AV  required  for  deflection  by 
several  times  the  value  calculated  using  a  two-body  approximation.  By  combining  a 
detailed  understanding  of  the  AV  requirements  for  NEO  deflection  with  the  analysis  of 
energy  coupling  by  Ahrens  and  Harris  (Ahrens,  1994,  p.  897),  an  estimate  of  the 
capability  of  current  impulsive  technologies  can  be  made.  It  is  important  to  remember 
that  the  analysis  presented  in  Chapter  VI  is  not  precise.  Many  approximations  and 
assumptions  were  combined  to  provide  an  educated  estimate  of  a  single  mission 
capability.  Significant  amounts  of  testing,  including  an  actual  NEO  deflection,  will  be 
required  to  build  upon  and  confirm  this  analysis.  Also,  there  were  several  assumptions 
made  in  the  patched-conic  analysis  for  simplification.  While  the  results  are  accurate,  and 
give  good  insight  into  the  dynamics  of  a  deflection  mission,  it  must  be  remembered  that 
the  numbers  are  presented  for  NEOs  in  elliptical  orbits  that  are  co-planar  with  the  Earth's 
orbit.  The  need  to  extend  the  analysis  into  the  third  dimension  is  discussed  as  a 
recommendation  for  further  study. 

B.  RECOMMENDATIONS  FOR  FURTHER  STUDY 


1.  Three-Body  Truth  Model 

A  significant  amount  of  work  has  now  been  accomplished  towards  rmderstanding 
the  orbital  mechanics  that  govern  a  NEO  deflection  mission.  The  best  way  to  test  this 
analysis,  barring  an  actual  NEO  deflection  test  mission,  would  be  to  build  a  computer 
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truth  model.  This  computer  code  could  model  the  three-body  dynamics  between  the  Sun, 
the  NEO,  and  the  Earth.  A  minimum  AV  from  the  above  analysis  can  be  applied  in  this 
model  to  confirm  that  the  NEO  is  diverted  as  intended. 


The  equations  of  motion  for  this  code  are  a  modified  version  of  the  solution  to  the 
classic  circular  restricted  three-body  problem.  In  the  classic  problem,  two  massive 
primaries  are  established  in  a  circular  orbit  about  their  center  of  mass  (Weisel,  1997,  p. 
278).  The  third  body  orbits  according  to  the  gravitational  influence  of  the  two  primaries. 
In  the  modified  circular  restricted  three-body  problem,  one  of  the  primaries  (the  Sim)  is 
significantly  more  massive  than  both  of  the  remaining  bodies.  As  such,  the  second 
primary  (the  Earth)  is  established  in  a  circular  orbit  not  around  the  center  of  mass  of  the 
primaries,  but  the  center  of  mass  of  the  Sun.  The  geometry  is  shown  in  Figure  7. 1 . 


The  ij  -coordinate  frame  rotates  about  the  Sun's  center  of  mass  at  a  rate  (D  equal  to 
the  rotation  rate  of  the  Earth  about  the  Sun.  The  Earth  lies  along  the  i  -axis  at  a  fixed 
distance  x, . 


First,  consider  the  equations  of  motion  as  governed  by  Newton's  Laws  of  Motion 
and  Gravity  (Wiesel,  1997,  p.  35) 

(7.1) 


m, 


pi  ^jt 


Applying  Equation  7.1  to  the  three  body  problem  shown  in  Figure  7.1  gives 


m^r2  = 


^3 


(7.2) 


The  zero  in  the  first  term  of  Equation  12  comes  fi'om  the  fact  that  the  Sun  is 
located  at  the  inertial  reference  point.  Dividing  through  by  mA  and  performing  some 
vector  algebra  gives 


1*2  = 


Gm^_ 

3  *2  3  *3 

h  h 


(7.3) 


The  choice  of  r2  rather  than  r3  for  the  equations  of  motion  now  becomes  apparent. 
By  using  r2,  the  mass  of  the  NEO  cancels  out  while  the  mass  of  the  two  primaries  (the 
Sun  and  the  Earth)  remains.  Thus  Equation  7.3  explains  the  gravitational  force  of  the  Sun 
and  Earth  applied  to  a  NEO  per  unit  mass. 

Some  manipulation  will  provide  the  acceleration  in  terms  of  the  coordinate  system 
components 

_  G^  j  2  ^  ^ i + y (7.4) 


m 


nia 


=  -G  ^^  +  ^  i-G  j 


^3 


(7.5) 


A  useful  result  would  be  to  reduce  Equation  7.5  to  two  variables  (X2  and  y2).  This 
can  be  accomplished  by  use  of  some  vector  relations. 


r,+r3=r2 


(7.6) 


r3=ra-r, 


(7.7) 


Thus 
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^3*  +  )'3  j  =  (^2  -  ^  1 )  >  +  ( 3^2  -  3’l  ) j 


(7.8) 


In  the  i  -component,  X2  is  a  variable  but  Xi  is  a  known  fixed  distance  of  1  AU. 
The  j -component  reduces  to  y2  because  yi  is  defined  to  be  equal  to  zero.  Placing  these 
results  into  Equation  7.5  gives 

(7.9) 

The  two  radii  can  be  expressed  in  terms  of  X2,  yi,  and  the  known  variable  Xj  as 

follows 

Tj  =-Jx2  +  y2  (7.10) 

r,=^{x,-lAvf  +  yl  (7.11) 

Equation  7.9  is  the  equation  of  motion  of  the  NEO  in  the  three-body  coordinate 
system.  Now  it  is  necessary  to  investigate  the  dynamics  of  the  system.  Beginning  with 
the  definition  of  the  r2  vector,  the  velocity  can  be  derived  as  follows: 

r2=^i  +  )’2j  (7.12) 

1-2  =  +  xji  +y2  j+  yii  (7.13) 

The  Transportation  Theorem  (Greenwood,  1988,  p.  49)  must  be  applied  because 
the  coordinate  system  is  rotating  about  the  k  -axis.  This  Theorem  states  that 


i  =  a*:xi=ttg  (7.14) 

and 

j=dkxj  =  -(oi  (7.15) 

Substitution  of  Equations  7.14  and  7.15  into  Equation  7.13  yields 

Tj  =i^i  +  X2(a)i)+y2j+3'2(-<»>)  (7-16) 

=  (i2-<wy2)l+(j2+«»^2)j  (7-17) 

The  velocity  can  now  be  differentiated  to  give  the  acceleration 

fj  ={x2  -a)y2)i+(i2-<«y2)i+(y2+aw:2)j+(y2  +  an:2)j  (7.18) 
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=  -o>j2)»+(^2  -«y2)(4)+(j2  +o«2)j+(>'2  +Cta2)(-a)i) 

fj  =[x2  “2(uy2  -(o^X2)i  +  {y2  +2(0X2  “®^^J2)j 


(7.19) 

(7.20) 


Equation  7.20  defines  the  dynamics  of  the  three-body  problem  in  the  rotating 
coordinate  system.  To  convert  to  a  state-space  form  useful  in  building  the  truth  model, 
the  components  of  Equations  7.9  and  7.20  are  equated.  This  process  yields  the  following 
final  relations 


With  Equations  7.21,  the  motion  of  the  NEO  has  been  described  using  two  second 
order  differential  equations.  Because  the  equations  are  non-linear,  the  best  solution  is  to 
convert  to  a  state-space  form  and  resolve  them  numerically  on  a  computer.  The  truth 
model  can  then  confirm  that  rs  equals  the  desired  minimum  as  calculated  in  the  EDAI 
code. 


2.  Three-Dimensional  Analysis 

The  assumption  of  co-planar  orbits  greatly  simplified  the  above  analysis,  but  it 
also  severely  limited  the  practicable  application  of  the  analysis.  Generalizing  the  theory 
to  three  dimensions  is  critical  in  application  to  the  majority  of  hazardous  NEOs. 

3.  Hyperbolic  Orbital  Analysis 

There  are  two  points  where  a  detailed  hyperbolic  analysis  can  be  added  to  the 
NEO  mitigation  analysis.  First,  many  NEOs  can  be  established  in  a  hyperbolic  orbit 
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about  the  Sun  before  impact  with  the  Earth.  These  orbits  were  not  included  in  this  thesis. 
Second,  the  effectiveness  of  an  impulse  applied  within  the  sphere  of  influence  of  the 
Earth  was  not  investigated  in  this  thesis.  Both  areas  are  important  for  a  thorough 
understanding  of  the  NEO  impact  problem. 

4.  EDAI  Code  Development 

The  EDAI  code  has  great  potential  to  be  very  useful  in  designing  a  mitigation  plan. 
As  currently  written,  the  code  generates  a  minimum  AV  required  to  generate  a  given  miss 
distance.  By  adding  the  analysis  included  in  Chapter  VI,  the  code  could  produce  a  more 
useful  result,  namely  the  ability  of  a  proposed  mission  to  effectively  mitigate  a  known 
threat.  In  a  fully  developed  code,  the  user  should  be  able  to  input  NEC's  orbital 
parameters  (including  the  inclination  with  respect  to  the  ecliptic),  the  NEC's  physical 
properties,  and  the  properties  of  proposed  mitigation  efforts.  The  code  would  then 
output  a  plot  of  the  mitigation  capability  against  the  diameter  or  mass  of  the  threat-NEC, 
as  in  Figure  6.1.  The  user  could  then  determine  if  the  mitigation  effort  is  sufficient  to 
handle  potential  errors  in  estimates  of  the  NEC's  physical  properties. 

As  the  first  user  not  involved  in  the  writing  of  the  code,  it  would  be  constructive 
to  provide  suggestions  for  improvement.  With  some  improvements  and  further  testing, 
the  code  could  be  v^ddely  used  in  researching  the  hazard  mitigation  problem. 

The  importance  of  an  accurate  initial  guess  was  discussed  in  Chapter  V.  As 
currently  written,  a  user  is  required  to  directly  modify  the  code  to  input  a  new  initial 
guess.  A  possible  improvement  would  be  to  allow  the  user  to  input  an  initial  guess. 
However,  this  would  require  a  fairly  well-educated  guess,  which  is  difficult  to  make 
without  a  thorough  understanding  of  the  problem.  Another  option  for  improvement 
would  be  to  build  a  reference  file  of  results  from  previous  successful  runs.  This  reference 
file  could  be  accessed  once  the  user  has  input  the  orbital  parameters  of  interest  to  get  an 
accurate  initial  guess. 
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Chapter  V  also  discusses  an  instability  induced  by  round-off  error.  The  solution 
to  this  problem  is  to  run  the  code  with  a  larger  minimum  separation  (e.g.,  10  R^)  and 
divide  the  result  to  achieve  the  desired  separation  (e.g.,  divide  by  10).  A  more  robust 
code  could  test  for  potential  round-off  instabilities  and  eliminate  this  complication  for  the 
user. 

It  would  be  helpful  to  have  the  ability  to  name  the  output  file  during  initialization 
while  the  code  is  actually  running.  This  would  eliminate  the  need  to  alter  the  code  to 
rename  the  output  file.  Also,  adding  the  NEO's  orbital  parameters  and  the  required 
separation  distance  to  the  output  file  could  help  to  avoid  confusion  between  different 
files. 

A  minor  bug  in  the  code  was  discovered  when  analyzing  NEO  orbits  with  a  semi¬ 
major  axis  equal  to  1  AU.  In  this  unique  case,  the  period  of  the  NEO  is  identical  to  the 
period  of  the  Earth.  Thus,  an  impact  would  occur  after  each  NEO  period.  The  EDAI 
code,  however,  only  calculates  the  minimum  AV  required  to  divert  the  NEO  for  a  specific 
impact  time.  The  code  does  not  recognize  that  impact  will  occur  after  each  orbital  period. 
Thus,  the  results  indicate  that  a  AV  applied  several  periods  before  impact  will  have  a 
greater  effect  than  a  AV  applied  within  an  orbital  period  of  impact.  In  reality,  the  results 
should  be  the  same  for  each  period. 

5.  Billiards  Scenario 

There  is  a  fourth  potential  impulse  concept  that  was  not  investigated  in  this 
analysis.  Melosh ,  Nemchinov,  and  Zetzer  (Melosh,  1994,  p.  1115)  suggest  diverting  a 
smaller  NEO  into  a  large  threatening  NEO.  It  could  be  useful  to  investigate  this  proposal 
using  techniques  developed  in  this  thesis. 
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6.  Verify  Results 


The  h5^othesis  discussed  in  Chapter  V  that  the  three-body  results  differ  from  the 
two-body  results  by  a  factor  exactly  equal  to  the  impact  radius  needs  to  be  verified. 
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APPElVfDIX  A.  V«>  PLOT  M-FILE 


%  Vasterth3.ni 

%  This  m-file  plots  the  magnitude  of  the  velocity  of  an  asteroid,  relative 
%  to  the  Earth,  at  the  impact  point.  The  variables  are  a_ast_sun  and  e_ast_sun. 
%  Assumptions  are: 

% 

%  -  The  Earth  is  in  a  circular  orbit 

clear 


%  First  set  some  constants 
mu_sun  =  1.32712438ell;  * 

mu_erth  =  3.986005e5;  % 

r_erth_sun  =  1.4959787e8;  % 

R_sun  =  6.96e5;  % 

a_min  =  (r_erth_sun  +  R_sun)/2;  % 
a_max  =  10*r_erth_sun;  % 

Re  =  6378;  % 

n_max  =  50; 


kmA3/s'^2 

kiiiA3/s'^2 

km 

km 

km 

km 

km 


%  The  first  FOR-loop  varies  e  from  0  to  1  at  a  step  of  0.1 
for  n  =  l:(n_max  -  1) 


%  The  next  (imbedded)  FOR-loop  varies  a  from  a_min  to  a_max 
for  step  =  l:(n_max  -  1) 
e(n,step)  =  n/n_max; 

a(n,step)  =  (step  -  l)*(a_max  -  a_min)/50  +  a_roin; 

%  Check  that  the  orbital  parameters  define  an  orbit  that  intersects  the  Earth's 
%  orbit. 

if  a(n,step)*(l-e(n,step))  >  r_erth_sun 

V_ast_erth(n,step)  =  0; 
bi(n,step)  =0; 

elseif  a(n,step)*(l+e(n,step))  <  r_erth_sun 

V_ast_erth(n,step)  =  0; 
bi(n,step)  =  0; 

else 

X  If  the  asteroid’s  orbit  will  intersect  the  Earth's  orbit,  compute  the  true  anomaly 
%  of  the  impact  and  the  flight  path  angle  at  impact. 

nu_imp_erth(n,step)  =  acos((a(n,step)*(l-e(n,step)>'2))/(e(n,step)*r_erth_sun)  -  1/ 
e(n,step)); 

garama_imp(n,step)  =  atan((e(n,step)*sin(nu_imp_erth(n,step)))/(l  +  e(n,step)* 
cos(nu_imp_erth(n , step)))) ; 
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%  Then  compute  the  magnitude  of  the  velocity  of  the  asteroid  with  respect  to  Earth. 

Vxl_ast_erth(n,step)  =  sqrt((2*mu_sun)/»’-e'’th_sun  -  mu_sun/a(n,step))* 
cos(gamma_impCn , step)) ; 

Vx_ast_erth(n,step)  =  Vxl_ast_erth(n,step)  -  sqrt(mu_sun/r_erth_sun); 
Vy_ast_erth(n,step)  =  sqrt((2*mu_sun)/r_erth_sun  -  mu_sun/a(n,step))* 
sin(gamma_imp(n , step)) ; 

V_ast_erth(n,step)  =  sqrt(Vx_ast_erth(n,step)'^2  +  Vy_ast_erth(n,step)A2); 
bi(n,step)  =  sqrt(C2*mu_erth)/CRe*V_ast_erth(n,step)'^2)  +1);  X  Re 

end 

[VmaxCn),i]  =  max(V_ast_erth(n, :)); 
aVmaxCn)  =  a(n,i)/r_erth_sun; 
eVmaxCn)  =  n/n_max; 

Cbimax(n),j]  =  niax(bi(n, :)); 
abimaxCn)  =  aCn, j)/r_erth_sun; 
ebimax(n)  =  n/n_max; 
end 


figure(l) 

col  ormapC  white') 
aplot  =  a/r_erth_sun; 
surf(aplot,e,V_ast_erth) 
grid  on 

xlabelC’a  (AU)') 
ylabel('e') 

zlabelC ' V_i_n_f_i_n_i_t_y  Ckm/s) ' ) 
titleC'V_i_n_f_i_n_i„t_y  vs.  a,e') 

figure(2) 

colormapC 'white ' ) 
surflCaplot,e,bi) 
grid  on 

xlabeK'a  (AU)’) 
ylabel(’e’) 

ZlabelC 'b_i  (Re)') 

titleC ’ Impact  Radius  (b_i)  vs.  a,e') 

Vi ew( -37. 5-125,  30) 

answer  =  input(' Generate  2-D  impact  radius  plot?  (y/n)  ','s’); 
while  answer=='y'; 

e_interest  =  input('At  which  eccentricity?  '); 

index  =  round(e_interest*n_raax); 

figure(3) 

plot(aplot(index, :) ,bi(index, :)) 
grid  on 

xlabel  ('a  (AU)') 
ylabel  (’bi  (Re)') 
e_title=num2str(e_interest) ; 
titleCC’bi  vs.  a  for  e  =  ',e_title]) 
print_ans  =  input(’Print  2-D  plot?  (y/n)  ','5'); 
if  print_ans  =«  'y' 
figure(3) 
print 
end 

answer  «  inputC Generate  2-D  impact  radius  plot?  (y/n)  ’,’s'); 
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APPENDIX  B.  MAPLE  ANALYSIS  FOR  V«,,  max 


File  name:  vbnd9.mws 

Last  Modified:  5  August  1997 

Written  by:  LT  Scott  Porter 

Purpose:  1 .  Determine  if  the  impact  velocity  between  the  Earth  and  an  asteroid  in  an  elliptical 

orbit  has  a  maximum. 

2.  If  a  maYimiim  exists,  find  a  fimction  of  an  asteoid's  semi-major  axis  and  it's 
eccentricity  which  results  in  the  maximum  impact  velocity  between  the  Earth  and  an  asteroid. 


This  analysis  uses  the  prinicples  of  orbital  mechanics,  vector  algebra,  and  siirqjle  calculus  to  determine 
if  the  relative  velocity  between  an  asteroid  and  the  Earth  at  impact  has  a  maximum.  The  asteroid  is 
assumed  to  be  in  an  elliptical  orbit  which  is  co-planar  with  the  Earth's  orbit.  The  Earth  is  assumed  to 
be  in  a  circular  orbit. 

>  restart:  with(linalg)  :  readlib. (isolate)  : 

Warning,  new  definition  for  norm 
.  Warning,  new  definition  for  trace 

Express  the  magnitude  of  the  velocity  of  the  Earth  and  the  asteroid. 

>  v[E] (a) :=sqrt (mu [sun] /r [E , sun] ) ; 

>  v[A,sun] (a) :=sqrt (2*mu[sun] /r [E,sun]  -  mu[sun]/a); 


Convert  these  velocities  to  vector  form.  The  coordinate  system  is  defined  with  the  X- Y  plane  in  the 
orbital  plane,  the  origin  at  the  center  of  mass  of  the  Earth,  and  the  Y-axis  continuosly  pointing  to  the 
center  of  mass  of  the  Sun.  Thus,  the  velocity  of  the  Earth  is  always  fixed  along  the  X-axis.  The 
velocity  of  the  asteroid  is  defined  by  the  flight  path  angle  (Gamma(a)),  the  angle  between  a  tangent  to  a 
circle  at  the  asteroid's  radius  and  the  velocity  vector. 

>  V[E ,sun]  (a)  :=vector (2 ,  [v[E]  (a)  , 0] )  ; 

>  V[A, sun] (a) :=vector (2 , [v [A, sun] (a) *cos (Gamma (a) ) , v [A, sun] (a) *sin (G 
amma  (a)  )  ] )  ; 
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Using  vector  algebra,  define  the  velocity  of  the  asteroid  with  respect  to  the  Earth  (V[A,E]).  V[A,E]  is 
defined  as  the  velocity  of  the  asteroid  with  respect  to  the  sun  (V[A,sun])  minus  the  velocity  of  the 
Earth  with  respect  to  the  sun  (V[E,sun]). 

>  V[A, E3  (a) : =evalm (V[A, sun] (a) -V[E , sun] (a) ) ; 

sin(r(fl)) 

Now  determine  the  orbital  position  and  flight  path  angle  of  the  asteroid  at  Earth  impact. 


>  nu  [I]  (a)  :  =arccos  (  (a*  (l-e''2)  /  (e*r  [E ,  sun] )  )  -  1/e); 

>  Gekinmal  (a)  :=arctan  (e*sin  (nu[I]  (a) )  /  (l+e*cos  (nu[I]  (a)  )  )  )  ; 


^  fl(l-e') 
v/ a)  :=  71  -  arccos  - 

I  ^^E.sun 


Substitute  the  impact  orbital  parameters  to  determine  the  velocity  of  the  asteroid  with  respect  to  the 
Earth  at  impact. 


>  V[A,E]  (a)  :  =subs  (Gamma  (a)  =Gammal  (a)  ,y  [A,E]  (a)  >  ; 
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I  '  / 

+  1 2  c7-  -  4  a  )  a  )  /  (a  { -r^_ +  2  a )  ), 


>  V3 [A,E]  (a)  :=vector (2 , [expand (V2 [A,E]  (a)  [1] )  , expand (V2 [A,E]  (a)  [2] )  j 

\  • 

t  f 

r 


"  >  V4[A,E]  (a)  :=vector(2,  [simplify (V3 [A, E]  (a)  [1])  , simplify (V3 [A, E]  (a)  [ 
2])])  ; 


Find  the  square  of  the  magnitude  of  the  velocity.  Using  the  square  of  the  magnitude,  rather  than  just 
the  magnitude,  simplifies  the  algebra  but  gives  the  same  final  result.  The  maximum  of  a  function  occurs 
at  the  same  place  as  the  maximum  of  it's  square. 
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>  Vniag2[A,El  (a)  :  =  (V3[A,Ej  (a)  [1]  ) ''2  +  V3[A,E]  (a)  [2]'^2 


Vmag2 1  a )  ;= 


^  ^  E.stm  ~  ^  jiin sj/n  ^  i^sim  ^  E.  sun  ^  \^sun  ’  E.  sun 

~  ^^sun  'E.sun  t^sun  ^E.sun  ~  ^  ^  ^~'’E.sun  +  2  fl)  ^£.ia„  ) 

Simplify  the  expression. 

>  Vmag[A,E] (a) :=simplify (") ; 

( 

2  - 

Vmag^  g{o)  '.—  —5  P5„„  ^  +  6  a  +  p^„„  /  e. ran 


+  2 


l^ran  ^  ^E.  sun  ~''  2  fl  ) 


a  r 


E.  sun 


-l^»,n  ( -1  +  ( -'e. ran  +  ^  ^ )  /  ( ^E. sun  «  ^ -'’£.5«n  +  2  a ) ) 

y 

combine ( ” ) ; 

^E.  ran  l^ran  ^  6  Pi„n  l^iun  ^  E.  jun  "*■  2  ( 

(f^ran  ^ ’"E.sun  “  ^  M-ran  ^  ^E. ran  ^  l^ran  ^  “  l^ran  ^E.ran  l^ran  ^  ^  ^E.ran~^l^iun  ^  ^  ) 

’’E.sun)''  /  ’^E.sun)/ (^E.««  «  ^ ^  fl  )  ) 

factor (") ; 

2  2^ 
sun  i^stm  ^  ^  M’sun  ^£:.  sun  ^  ^  ^£,  sun 

+  2  a  7 -p, J  (e  -  1 )  (e  +  1 )  +  2  a)"  )  /  (r^.„y  a  (-r£,,„„  +  2  a)) 

collect  ("  ,a)  ;  _ _ _ . 

^^«n  ’’E.sun  +  (“^  l^ran  ''E.ran*  +  2  7-^^^</n'  «  («  '  1  )  (^  +  1  )  (“^E.^nn  +  2  a)  r£.„„  )  « 
r  /  ? 


sun  e.sun  ''  «  sun  c..sun  y  » 

^V-sun’’E,stm)  /  (’'E.s.m  ^  (-’"E.sun  +  ^  ^)) 
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>  Viiiagsqrd[A,E]  (a)  :=sort(''}  ; 

Vmagscjni^  li a) -={6  a  /*£  >'e. sun 

a)/  ((2a-rt-^„)a 

’'E.sim  f 

Differentiate  the  function  for  the  magnitude  (squared)  and  set  the  result  equal  to  zero  to  find  the 
maximum. 


>  Vmaxl[A,E]  (a)  :  =siiaplify  (dif f  (Vmagsqrd [A, E]  (a)  ,a)-0)  ; 

Vmaxl  f  i.i a)  :=  -  (-2  |ij„„  o'  +  2  e~  +  a~  ''e.suu  ~  Ksfm  ^'E.stm  ^ 

-  J„„  '\J ~(e  —  \  )  {e  +  \  )  {2  a  -  lti„n  ^'e.suu  )  l^?Hn  ^  E. sun 

^  \  +  I  )  (  2  ■“  /  )  n  /  E.  sun  ^  ~  ^ 

>  Vtest[A,E]  (a)  :  =siiBplify  (diff  (Vmaxl  [A,E]  (a)  ,a)  )  ; 


I  Solve  this  function  for  the  semi-major  axis.  This  result  will  determine,  for  a  given  eccentricity  (e),  the 
i  semi-major  axis  which  will  generate  the  maximum  relative  velocity  between  the  asteroid  and  the  Earth 

at  impact. 


> 


isolate (Vmaxl [A, E]  (a) ,a)  ; 


,2  1/3 

(-(-I  +e‘)  )  ^E.sun 


I  >  subs (a=amax, ") ; 

! 


,  2  1/3 

^  E.  sun 

amax  =  ; 

-1  +e‘' 


r  >  assign (”) ; 


To  further  simplify  the  relation,  use  r[E,sun]=l  AU, 

>  aAU:  =s\ibs  (r  [E  ,  sun]  =1 ,  amax)  ; 
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aAU-= 


,2  (;3 

(-(-1  +e-)) 

-1  +  e- 


>  Vtest:2tA,E]  (a)  :  =subs  (a=aAU,  Vtest  (A,  E]  (a))  ; 


1 

Vtest2  ,  £-(«):=--  2  —  7  -  2  — 

2  -1  +  e"  +e 


(-(-1  +£?■))  i^siin''E.su 


(-1 


i-l+e-) 


(  2  1/3  ^2 

(-{-1  +£?-)  )  ^  ^  ,  .  2  1/3  : 

(S—  \)i6+\)2  ,  ~  ^E.sun  ^  (  1  +  )  l^sun  ^E.sun 

_ \ _ -1  +  g~ _ / _ . 

-1  +t;- 


(-l+e-)u. 


(-(-1+g-))  ,  2%'^^  2 

(e-l)(e  +  l)  2  -^Esun  (-(-1+c)) 

V  -1  +  e“  J 


-1  +e“ 


>  simplify  (")  ; 


1  4  '*£.  sun  ^  f*  -  2  -  4  r^_  +  (  -(  - 1  +  t?'  )  )  /  £.  ) 


e‘)  (-(-l+e‘)  ) 


i-l+e-) 


>  Slabs  (r  [E ,  sun]  =1 ,  " )  ; 


1  (4  +  2  U^„  e~-2  u  -  4  +  (-(-1  +  c" )  ) 


(-2(-(-l +e-)  )  -1 +e“)  (-(-1 +e‘)  ) 


>  simplify (">  ; 


2  7^ 


%1  :=-• 


,21/3  ,2  ,21/3  2 

(-2(-(-l +t?‘)  )  -  1 +^‘)  (-(-1 +e')  ) 


(-1  +e- 


>  eval (subs (mu [sun] =1 , ” ) ) ; 

1  4V^e'  +  2g^-2-4V^+(-(-l  +e-)~)~'' 

2  7^ 


=  0 


2  1/3 


%1  :=- 

>  simplify (") ; 


(-2(-(-l +g-)~)  '  -l+g-)~(-(-l +g~)  ) 

2 

(-i+g') 


2  1/2 


1  47^g~  +  2g~-2-47%l  +(-(-l+g~)~) 


=  0 


2  7^ 

,21/3  ,2  ,2  1/3 

{-2(-(-l+g“)  )  -1+g-)  (-(-l+g-)  ) 

(-1  +g‘) 


%1  :=-■ 


End  of  analysis. 
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APPENDIX  C.  INITIAL  GUESSES  FOR  EDAI  CODE 


As  with  any  numerical  optimization  algorithm,  the  EDAI  code  requires  an 
accurate  initial  guess  of  the  results.  An  inaccurate  guess  will  result  in  an  instability.  For 
this  analysis,  the  code  was  first  run  with  a  known  solution  used  to  provide  an  initial 
guess.  Then  the  NEO's  orbital  parameters  were  altered  slightly,  such  that  the  initial  guess 
from  the  previous  run  was  still  close  enough  to  the  correct  answer  to  provide  a  stable  run. 
This  process  was  continued  until  a  database  of  initial  guesses  was  built.  The  initial 
guesses  are  provided  here. 

Canonical  units  are  used  for  the  AV  and  tf  guesses.  The  units  are  defined  as 
follows: 


DU  =  Distance  Unit  =  Astronomical  Unit  (AU)  =  1 .495978e8  km 

TU  =  Time  Unit  =  —  year 
2nr 


e 

a(AU) 

Minimum 
Separation  (R^ 

AV||  (DU/TU) 

I  AVx(DU/rU) 

tf(TU) 

0.1 

1.0  i 

10 

1.062e-7 

I  4.142e-4  i 

2.843e-4 

0.2 

1.0 

10 

4.963e-5 

1  4.157e-4 

1.336e-4 

0.3 

1.0 

10 

-1.581e-5 

j  3.999e-4 

8.856e-5 

0.4 

1.0 

10 

-8.548e-5 

I  3.575e-4 

6.91 Oe-5 

0.5 

1.0 

10 

-1.4417e-4 

I  2.769e-4  | 

1  5.543e-5 

0.6 

1 

1 . . . ™50 

-7.926e-4 

1  8.303e-4  | 

1  1.674e-4 

r  r.6  1 

50 

-6.057e-4 

1  4.340e-4  i 

1  3.446e-7 

1  1-0  ^ 

1 

-4.221e-4 

3.440e-4 

1  -1.331e-4 

0.9 

1  1-0 

1  100 

-5.087e-4 

1  7.21 8e-4 

1  -4.042e-4 
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